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ABSTRACT 


The  application  of  coupled  Lagrangian-Eulerian  and  multimaterial  Eulerian  finite  element 
analysis  to  problems  of  interest  in  underwater  shock  research  was  investigated.  Analyses  were 
conducted  for  the  classical  problems  of  a  spherical  shell  and  an  infinite  cylinder  loaded  by  a 
plane  acoustic  step  wave,  and  for  the  expansion  and  collapse  of  explosion  generated  gas  bubbles 
in  deep  water.  Results  for  the  elastic  response  of  a  spherical  shell  and  an  infinite  cylinder  were 
found  to  be  in  excellent  agreement  with  the  analytical  solution  results,  and  results  for  the 
expansion  and  collapse  of  a  deep  explosion  generated  gas  bubble  in  the  absence  of  nearby 
boundary  surfaces  were  found  to  agree  well  with  experimental  results.  Several  analyses  were 
conducted  with  explosive  charges  detonated  at  various  distances  from  a  rigid  wall;  results  from 
these  analyses  were  in  qualitative  agreement  with  what  is  known  about  this  case,  and  serve  to 
characterize  the  behavior  of  the  resulting  bubble  in  terms  of  the  distance  to  the  wall. 


Aeaeaaion  Fas? 

'litis  Q’iiAkl 

TAIJ 

□ 

Unsimovisicad 

□ 

Jus  t  i  f  1  c  a  t  i  on _ 

By - - 

Distr 

liva-ti  aad'/oT 


TABLE  OF  CONTENTS 


I.  INTRODUCTION .  1 

II.  NUMERICAL  COMPUTER  CODE .  3 

III.  COUPLED  ANALYSIS  OF  CLASSICAL  PROBLEMS  .  7 

A.  SPHERICAL  SHELL/PLANE  STEP  WAVE .  7 

B.  INFINITE  CYLINDER/PLANE  STEP  WAVE .  13 

IV.  MULTIMATERIAL  EULERIAN  ANALYSIS  OF  BUBBLES .  21 

A.  DEEP  SPHERICAL  BUBBLE .  21 

B.  BUBBLES  NEAR  RIGID  BOUNDARIES .  29 

V.  CONCLUSION  . 53 

APPENDIX  A:  SPHERICAL  SHELL/PLANE  STEP  WAVE  INPUT  ....  55 

APPENDIX  B:  INFINITE  CYLINDER/PLANE  STEP  WAVE  INPUT  ...  59 

APPENDIX  C;  A  POLYNOMIAL  STATE  EQUATION  FOR 

SEAWATER .  61 

APPENDIX  D;  DEEP  SPHERICAL  BUBBLE  INPUT .  65 

APPENDIX  E:  BUBBLES  NEAR  RIGID  BOUNDARIES-SHAPES  ....  69 

APPENDIX  F:  BUBBLES  NEAR  RIGID  BOUNDARIES-INPUT .  97 

LIST  OF  REFERENCES  .  115 

INITIAL  DISTRIBUTION  LIST . 118 


V 


LIST  OF  FIGURES 


Figure  1 .  Eulerian  Element  Control  Volume  Changes  Due  to  Presence 

of  Coupling  Surface  With  "General  Coupling"  Method .  4 

Figure  2.  Blending  of  Eulerian  Elements  to  Prevent  Small  Element  From 
Controlling  Stable  Time  Step  Size  When  Using  "General 

Coupling"  Method .  5 

Figure  3.  Comparison  Between  Eulerian  Mesh  With  "General  Coupling" 

and  "ALE  Coupling”  Methods  .  6 

Figure  4.  Spherical  Shell/Plane  Step  Wave  Problem  Geometry .  8 

Figure  5.  Steel  Shell  Elements,  Dummy  Elements,  and  Resulting  Closed 
Coupling  Surface  for  Spherical  Shell/Plane  Step  Wave 

Problem .  9 

Figure  6.  Fluid  Mesh  for  Spherical  Shell/Plane  Step  Wave  Problem  ...  11 

Figure  7.  Size  and  Position  of  Spherical  Shell  Relative  to  Eulerian  Fluid 

for  Spherical  Shell/Plane  Step  Wave  Problem  .  12 

Figure  8.  Non-Dimensional  Incident  Pressure  Independent  Radial 
Velocity  Versus  Non-Dimensional  Time  for  Spherical 

Shell/Plane  Step  Wave  Problem .  14 

Figure  9.  Infinite  Cylinder/Plane  Step  Wave  Problem  Geometry  .  15 

Figure  10.  Structural  and  Dummy  Elements  in  Finite  Element  Model  for 

Infinite  Cylinder/Plane  Step  Wave  Problem .  16 

Figure  1 1.  Fluid  Mesh  for  Infinite  Cylinder/Plane  Step  Wave  Problem  .  .  18 

Figure  12.  Size  and  Position  of  Infinite  Cylinder  Relative  to  Eulerian 

Fluid  for  Infinite  Cylinder/Plane  Step  Wave  Problem .  19 

vii 


Figure  13.  Non-Dimensional  Incident  Pressure  Independent  Radial 
Velocity  Versus  Non-Dimensional  Time  for  Infinite 

Cylinder/Plane  Step  Wave  Problem . 

Figure  14.  Deep  Spherical  Bubble  Problem  Geometry  . 

Figure  15.  Finite  Element  Model  for  Deep  Spherical  Bubble  Problem  .  . 

Figure  16.  Radius  Versus  Time  Behavior  for  Deep  Spherical  Bubble 

Problem  (0.2990  kg  TNT  Charge  at  178.6  m  Depth  in 

Seawater) . 

Figure  17.  Geometry  of  Bubble  Near  Rigid  Boundary  Problems  . 

Figure  1 8.  Overall  Model  Geometry  for  Quarter  Symmetry  Axisymmetric 

Finite  Element  Model  For  Free-Field  Analysis . 

Figure  19.  Quarter  Symmetry  Axisymmetric  Finite  Element  Model  for 

Free-Field  Analysis . 

Figure  20.  Close  Up  View  of  Finite  Element  Model  in  Area  Near  Charge 

for  Free-Field  Analysis  . 

Figure  2 1 .  Finite  Element  Model  Geometry  for  Rigid  Wall  at  1 .000  Free- 

Field  Radii  From  Center  of  Charge . 

Figure  22.  Close  Up  View  of  Finite  Element  Model  in  Area  Near  Charge 
for  Rigid  Wall  at  1.0(K)  Free-Field  Radii  From  Center  of 

Charge  . 

Figure  23.  Finite  Element  Model  Geometry  for  Rigid  Wall  at  4.015  Free- 

Field  Radii  From  Center  of  Charge . 

Figure  24.  Finite  Element  Model  for  Rigid  Wall  at  4.01 5  Free-Field  Radii 

From  Center  of  Charge . 

Figure  25.  Close  up  View  of  Finite  Element  Model  in  Area  Near  Charge 
for  Rigid  Wall  at  4.015  Free-Field  Radii  From  Center  of 
Charge  . 


20 

22 

26 

27 

30 

32 

33 

35 

36 

37 

39 

40 

41 


viii 


Figure  26.  Equivalent  Spherical  Radius  vs  Time  for  Explosion  Gas 
Bubbles  Near  a  Rigid  Wall  . 

Figure  27.  Non-Dimensional  First  Oscillation  Period  vs  (Non-Dimensional 
Standoff  Distance  From  Rigid  Wall) ' . 

Figure  28.  Displacement  of  Center  of  Mass  vs  Time  for  Explosion  Gas 
Bubbles  Near  a  Rigid  Wall  . 

Figure  29.  Maximum  and  Minimum  Equivalent  Spherical  Radii  Versus 
(Non-Dimensional  Standoff  Distance  From  Rigid  Wall)  ’  .  .  . 

Figure  30.  Displacement  of  Bubble  Center  of  Mass  at  First  Minimum 
and  Peak  Bubble  Velocity  Versus  (Non-Dimensional  Standoff 
Distance)  '  . 

Figure  31.  Bubble  Shape  for  Rigid  Wall  at  h*=4.015  at  Time  t=0.50T  .  . 

Figure  32.  Bubble  Shape  for  Rigid  Wall  at  h*=4.015  at  Time  t=0.80T  .  . 

Figure  33.  Bubble  Shape  for  Rigid  Wall  at  h*=4.015  at  Time  t=0.90T  .  . 

Figure  34.  Bubble  Shape  for  Rigid  Wall  at  h*=4.015  at  Time  t=0.95T  .  . 

Figure  35.  Bubble  Shape  for  Rigid  Wall  at  h*=4.()15  at  Time  t=0.98T  .  . 

Figure  36.  Bubble  Shape  for  Rigid  Wall  at  h*=4.015  at  Time  t=0.99T  .  . 

Figure  37.  Bubble  Shape  for  Rigid  Wall  at  h*=4.015  at  Time  t=1.00T  .  . 

Figure  38.  Bubble  Shape  for  Rigid  Wall  at  h*=2.008  at  Time  t=0.50T  .  . 

Figure  39.  Bubble  Shape  for  Rigid  Wall  at  h*=2.008  at  Time  t=0.80T  .  . 

Figure  40.  Bubble  Shape  for  Rigid  Wall  at  h*=2.008  at  Time  t=0.90T  .  . 

Figure  41.  Bubble  Shape  for  Rigid  Wall  at  h*=2.008  at  Time  t=0.95T  .  . 

Figure  42.  Bubble  Shape  for  Rigid  Wall  at  h*=2.()08  at  Time  t=0.98T  .  . 

Figure  43.  Bubble  Shape  for  Rigid  Wall  at  h*=2.008  at  Time  t=0.99T  .  . 

Figure  44.  Bubble  Shape  for  Rigid  Wall  at  h*=2.008  at  Time  t=1.00T  .  . 

Figure  45.  Bubble  Shape  for  Rigid  Wall  at  h*=1.374  at  Time  t=0.50T  .  . 

Figure  46.  Bubble  Shape  for  Rigid  Wall  at  h*=1.374  at  Time  t=0.80T  .  . 


43 

44 

46 

50 


51 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 


IX 


Figure  47.  Bubble  Shape  for  Rigid 
Figure  48.  Bubble  Shape  for  Rigid 
Figure  49.  Bubble  Shape  for  Rigid 
Figure  50.  Bubble  Shape  for  Rigid 
Figure  51 .  Bubble  Shape  for  Rigid 
Figure  52.  Bubble  Shape  for  Rigid 
Figure  53.  Bubble  Shape  for  Rigid 
Figure  54.  Bubble  Shape  for  Rigid 
Figure  55.  Bubble  Shape  for  Rigid 
Figure  56.  Bubble  Shape  for  Rigid 
Figure  57.  Bubble  Shape  for  Rigid 
Figure  58.  Bubble  Shape  for  Rigid 


Wall  at  h*=l  .374  at  Time  t=0.90T  . 
Wall  at  h*=1.374  at  Time  t=().95T  . 
Wall  at  h*=1.374  at  Time  t=0.98T  . 
Wall  at  h*=1.374  at  Time  t=0.99T  . 
Wall  at  h*=1.374  at  Time  t=1.00T  . 
Wall  at  h*=1.00()  at  Time  t=0.50T  . 
Wall  at  h*=1.000  at  Time  t=0.80T  . 
Wall  at  h*=1.00()  at  Time  t=0.90T  . 
Wall  at  h*=1.000  at  Time  t=0.95T  . 
Wall  at  h*=1.0(K)  at  Time  t=0.98T  . 
Wall  at  h*=1.0{X)  at  Time  t=0.99T  . 
Wall  at  h*=1.0(X)  at  Time  t=:1.00T  . 


85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 


X 


I.  INTRODUCTION 


This  report  describes  work  ongoing  at  the  Naval  Postgraduate  School  to 
apply  coupled  Lagragian-Eulerian  and  multimaterial  Eulerian  finite  element  analysis 
techniques  to  problems  in  underwater  shock  research.  Previous  underwater  shock 
research  at  the  Naval  Postgraduate  School  has  concentrated  on  experimentation,  e.g. 
by  Jones  and  Shin  [Ref.  1],  and  doubly  asymptotic  approximation  boundary 
element  techniques,  sometimes  combined  with  experimentation,  e.g.  by  Fox  et.  al 
[Ref.  2],  Nelson  et.  al.  [Ref.  3],  Kwon  et.  al  [Ref.  4],  and  Chisum  [Ref.  5]. 

However,  for  certain  types  of  problems  doubly  asymptotic  approximation 
boundary  element  techniques  have  not  been  advanced  to  the  point  that  they  can 
provide  useful  results.  This  is  true  in  particular  for  problems  involving  the 
oscillation  of  gas  bubbles  generated  by  underwater  explosions.  This  area  is  of 
considerable  interest,  as  the  pressure  pulsations  produced  by  such  bubbles  can 
under  certain  circumstances  produce  significant  whipping  of  nearby  marine  or 
submarine  structures. 

In  addition,  there  are  certain  practical  constraints  on  the  ability  to  conduct 
experimentation  to  determine  structural  responses.  Full  scale  experimentation  is 
extremely  expensive,  and  certain  physical  phenomena  related  to  explosion  gas 
bubbles  cannot  be  scaled  in  a  practical  experimental  setup. 

In  light  of  these  factors,  and  taking  into  consideration  the  ongoing  advances 
in  computer  capabilities  and  the  recent  availability  of  advanced  finite  element 
programs  which  can  efficiently  calculate  the  fluid-structure  interaction  between 
Eulerian  and  Lagrangian  materials  and  are  capable  of  dealing  with  several  different 
Eulerian  materials  in  the  same  problem,  a  more  basic  approach  might  allow 
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solution  of  heretofore  unsolvable  problems.  In  this  approach,  each  of  the  materials 
in  an  underwater  shock  problem  are  modeled  in  the  most  advantageous  way  for  that 
class  of  material;  fluid  media  and  explosives  using  Eulerian  elements,  and  structural 
materials  using  Lagrangian  elements. 

One  advantage  of  this  approach  is  that  there  are  few  approximations 
involved;  the  resulting  solution  can  essentially  be  made  as  accurate  as  the 
discretization  allowed  by  the  available  computational  resources  and  the  certainty 
with  which  the  properties  of  the  materials  involved  are  known  will  permit. 

This  approach  also  overcomes  the  problems  involved  with  modeling  all  of 
the  material  in  a  problem  with  Lagrangian  elements,  which  in  an  underwater  shock 
problem  quickly  become  so  distorted  that  the  stable  time  step  size  approaches  zero 
and  the  time  to  compute  a  solution  out  to  near  steady  state  approaches  infinity. 
Nor  is  the  approach  of  modeling  all  of  the  material  in  an  underwater  shock  problem 
using  Eulerian  materials  generally  practical,  as  this  approach  requires  that  an 
extremely  large  number  of  Eulerian  elements  be  used  in  order  to  accurately  capture 
the  response  of  structural  materials  in  the  problem,  which  is  usually  the  primary 
item  of  interest. 

By  using  a  finite  element  code  which  contains  both  Lagrangian  and  Eulerian 
processors  and  a  method  for  computing  the  fluid-structure  interaction  at  the 
interface  between  Lagrangian  and  Eulerian  materials,  the  advantages  of  both  types 
of  analysis  are  realized  and  the  shortcomings  associated  with  attempting  to  use  one 
or  the  other  alone  are  eliminated.  This  is  the  approach  we  have  chosen  to  pursue. 
Initial  results  from  this  approach  have  been  very  promising  [Ref.  6]. 
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II.  NUMERICAL  COMPUTER  CODE 


The  finite  element  program  used  for  the  results  described  in  this  report  is 
MSC/DYTRAN  [Ref.  7],  a  three  dimensional  finite  element  program  available  from 
the  MacNeal-Schwendler  Corporation.  This  program  was  developed  by  combining 
and  extending  two  other  computer  programs,  MSC/DYNA  [Ref.  8]  and 
MSC/PISCES  [Ref.  9].  Both  of  these  programs  have  a  proven  record  for  the 

analysis  areas  for  which  they  were  developed. 

Like  MSC/DYNA,  MSC/DYTRAN  is  capable  of  handling  non-linear,  large 
strain  structural  response  problems.  MSC/DYTRAN  is  also  capable  of  solving 
problems  involving  Lagrangian-Lagrangian  two  surface  (contact-impact)  and  single 
surface  (folding)  problems.  A  complete  constitutive  model  can  be  defined  in  terms 
of  an  equation  of  state,  a  shear  model,  a  yield  model,  a  failure  model,  and  a  spall 
model. 

Multimaterial  Eulerian  processors  in  MSC/DYTRAN  allow  for  up  to  9 
different  eulerian  materials  to  be  present  in  a  given  problem.  In  addition,  two 
different  methods  are  available  to  provide  for  calculation  of  the  fluid- structure 
interaction  between  Lagrangian  and  Eulerian  materials. 

In  the  "General  Coupling"  method,  the  Lagrangian  and  Eulerian  meshes  are 
geometrically  independent,  and  interact  via  a  coupling  surface  attached  to  the 
Lagrangian  structure.  This  method  requires  that  the  coupling  surface  form  a  closed, 
simply  connected  volume,  on  one  side  (inside  or  outside)  of  which  the  Eulerian 
elements  are  "void"  (contain  no  material).  The  deformable  coupling  surface  cuts 
across"  Eulerian  elements,  changing  their  control  volume  and  surface  areas.  This 
is  illustrated  in  Figure  1,  for  two  dimensional  Eulerian  elements.  To  prevent  the 
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Coupling 

Surface 

Uncovered  Side  /  Covered  Side 

Eulerian 
Elements 


Figure  1.  Eulerian  Element  Control  Volume  Changes  Due  to  Presence  of 
Coupling  Surface  With  "General  Coupling"  Method 

stable  time  step  size  from  being  controlled  by  very  small  Eulerian  control  volumes 
formed  by  the  coupling  surface,  elements  for  which  the  ratio  between  the  covered 
(void)  volume  fraction  and  the  initial  volume  is  less  than  a  user  modifiable  blend 
parameter  are  combined  with  adjacent  elements  to  form  larger  elements;  this  is 

illustrated  in  Figure  2. 

The  other  method  provided  by  MSC/DYTRAN  for  coupling  of  Lagrangian 
and  Eulerian  materials  is  "Arbitrary  Lagrange-Euler"  (ALE)  coupling.  In  this 
method,  the  fluid  and  structural  mesh  geometries  are  not  independent.  Instead,  the 
interface  surface  between  the  Lagrangian  and  Eulerian  elements  is  actually 
composed  of  the  union  of  the  faces  of  these  eleinents.  As  this  interface  is 
deformed  during  deformation  of  the  Lagrangian  structure,  Eulenan  gnd  points 
which  are  attached  to  this  also  move.  To  keep  the  geometry  of  the  Eulerian  mesh 
relatively  "nice",  other  Eulerian  grid  points  away  from  the  coupling  surface  can  be 
allowed  to  move,  e.g.  towards  the  center  of  their  nearest  neighbors.  Note  that  in 
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/ 

Figure  2.  Blending  of  Eulerian  Elements  to  Prevent  Small  Element  From 
Controlling  Stable  Time  Step  Size  When  Using  "General  Coupling"  Method 

this  method,  the  Eulerian  mesh  is  not  stationary.  However,  the  motion  of  the 
Eulerian  mesh  is  purely  geometrical;  the  velocity  of  material  through  this  mesh  is 

independent  of  the  motion  of  the  mesh. 

Figure  3  illustrates  the  difference  between  these  two  fluid-structure  coupling 

methods,  in  two  dimensions. 
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General  Coupling 

(Eiilerian  mesh  doesn’t  change) 


t  =  0 


7 

7 

7 

□ 

Coupling 

Surface 


t  =  1 


Surface 


ALE  Coupling 

(Eulerian  mesh  deforms  "vrith  structure) 


Figure  3.  Comparison  Between  Eulerian  Mesh  With  "General  Coupling"  and 
"ALE  Coupling"  Methods 
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III.  COUPLED  ANALYSIS  OF  CLASSICAL  PROBLEMS 


In  order  to  examine  the  performance  of  using  coupled  Lagrangian-Eulerian 
finite  element  analysis  for  underwater  shock  problems,  two  classical  problems  for 
which  analytical  solution  data  is  available  were  analyzed.  These  analyses  examined 
the  elastic  response  of  a  spherical  shell  and  an  infinite  cylinder  to  loading  from  a 
plane  step  wave  propagating  through  an  acoustic  media. 

Huang  has  solved  these  problems  analytically,  using  a  direct  inverse  Laplace 
transform  of  a  finite  number  of  terms  of  the  infinite  series  expansion  of  the 
equations  for  the  respective  shells  [Refs.  10  and  11].  For  our  finite  element 
analyses,  the  same  material  properties,  parameters,  and  non-dimensionalization 
procedures  used  by  Huang  in  his  analyses  were  used. 

A.  SPHERICAL  SHELL/PLANE  STEP  WAVE 

Figure  4  shows  the  geometry  of  the  spherical  shell/plane  step  wave 
problem.  The  material  properties  and  parameters  used  for  this  problem  were; 


Shell  Material; 

Young’s  Modulus  for  Steel;  30  x  10  psi 

Poisson’s  Ratio  for  Steel;  0.3 

Density  of  Steel;  Ibm/ft 

Shell  Thickness  to  Radius  Radio;  0.02 

Fluid;  Water  ^ 

Water  Density;  62.4  Ibm/ft 

Water  Acoustic  Wave  Speed;  4794  ft/s 


The  problem  was  non-dimensionalized  using  the  radius  of  the  sphere  as  the 
characteristic  length,  the  time  for  an  acoustic  wave  to  transit  one  radius  as  the 
characteristic  time,  and  the  bulk  modulus  of  water  as  the  characteristic  pressure. 
A  bulk  modulus  equation  of  state  was  used  to  model  the  water  in  this  problem,  and 
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Figure  4.  Spherical  Shell/Plane  Step  Wave  Problem  Geometry 

a  small  incident  pressure  wave  magnitude  (1  x  10'^  bulk  modulus  units)  was  used 
to  keep  deformations  small  enough  for  the  elastic  assumption  to  be  valid. 

For  our  finite  element  model,  a  quarter  symmetry  model  was  used.  An 
elastic  material  model  consisting  of  150  quadrilateral  Lagrangian  shell  elements 
was  used  to  model  one  quarter  of  the  spherical  steel  shell.  A  single  constraint  set 
was  used  to  constrain  the  appropriate  translational  and  rotational  degrees  of 
freedom  of  grid  points  lying  on  symmetry  planes. 

MSC/DYTRAN’s  "General  Coupling"  fluid-structure  interaction  method, 
in  which  the  Lagrangian  and  Eulerian  meshes  are  independent  and  interact  via  a 
coupling  surface,  was  used  for  this  problem.  This  method  requires  that  the 
coupling  surface  form  a  closed,  simply  connected  volume;  for  simplicity,  this 
closed  volume  was  generated  by  using  450  dummy  elements  in  addition  to  the  150 
Lagrangian  shell  elements  used  to  model  the  steel  shell.  The  Lagrangian  (steel) 
shell  elements,  the  dummy  elements,  and  the  resulting  closed  coupling  surface  are 
illustrated  in  Figure  5. 
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steel  Shell  Elements 
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Figure  5.  Steel  Shell  Elements,  Dummy  Elements,  and  Resulting  Closed 
Coupling  Surface  for  Spherical  Shell/Plane  Step  Wave  Problem 
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Because  only  a  finite  volume  of  fluid  material  can  be  modeled  using  this 
approach,  it  was  decided  to  construct  a  model  for  which  the  solution  would  be 
unaffected  by  reflection  from  the  boundaries  of  the  fluid  volume  for  times  less  than 
six  radius  transit  times.  The  block  of  water  modeled  is  thus  a  rectangle  bounded 
by  the  planes  x=0  and  x=4,  y=0  and  y=4,  and  z=4  and  z=-4,  where  the  point 
(0,0,0)  represents  the  center  of  the  sphere  and  units  are  in  terms  of  the  radius  of 
the  sphere.  Every  point  on  the  shell  is  thus  at  least  three  radii  away  from  a 
boundary,  and  since  acoustic  waves  travel  one  shell  radius  per  radius  transit  time, 
no  boundary  reflection  reaches  the  shell  for  six  radius  transit  times.  The  fluid 
mesh  used  consists  of  65,536  cubical  Eulerian  elements;  the  length  of  each  side  of 
each  element  is  1/8  radii.  Figure  6  shows  the  fluid  mesh  used  for  this  problem. 

All  boundaries  of  this  fluid  volume  were  left  with  a  default  wall  boundary 
condition  except  the  boundary  at  z=4  radii;  this  boundary  was  given  a  "flow- 
boundary  condition,  with  a  pressure  of  0.001  bulk  moduli!  and  a  particle  velocity, 
determined  from  the  one  dimensional  wave  equation 

p=pcu 

of  0.001  times  the  acoustic  wave  speed  in  water,  in  the  -z  direction.  Initial 
conditions  were  imposed  on  all  of  the  Eulerian  elements  such  that  all  elements 
between  the  z=4  radii  and  z=l  radii  planes  had  an  initial  pressure  of  0.001  bulk 
moduli!  and  a  particle  velocity  of  0.001  times  the  acoustic  wave  speed  in  water  in 
the  -z  direction,  and  all  elements  between  the  z=l  and  z=-4  radii  planes  had  zero 
initial  pressure  and  particle  velocity.  Time  t=0  for  the  finite  element  analysis  thus 
corresponds  to  the  instant  when  the  plane  step  wave  first  touches  the  sphere,  at  the 

point  (0,0,1). 

The  relationship  between  the  size  of  the  spherical  shell  and  the  fluid  volume 
modeled  in  this  problem  is  illustrated  in  Figure  7,  from  two  different  viewpoints. 
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Mesh  for  Spherical  Shell/Plane  Step  Wave  Problem 
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For  clarity,  only  the  outline  of  the  fluid  block  that  was  modeled  is  shown. 

The  steel  shell  element,  dummy  element,  and  fluid  element  geometry  for  this 
problem  were  all  created  using  the  finite  element  preprocessing  program  LS- 
INGRID  [Ref.  12].  A  small  FORTRAN  translating  program  was  then  used  to 
convert  the  output  files  from  LS-INGRID  to  NASTRAN  compatible  format  files, 
which  were  then  read  into  MSC/XL  [Ref.  13].  Symmetry  constraints,  the  "flow 
boundary  condition,  and  material  properties  were  then  added  using  MSC/XL,  and 
a  single  "bulk  data"  file  was  written  out.  The  LS-INGRID  input  files,  FORTRAN 
geometry  translating  program,  and  MSC/XL  input  data  stream  used  to  create  the 
model  for  this  problem,  together  with  the  MSC/DYTR  AN  input  file  used  to  control 
the  analysis,  are  provided  in  Appendix  A.  The  information  provided  in  Appendix 
A  is  sufficient  to  duplicate  our  analysis. 

The  resulting  transient  solution  for  the  radial  velocity  of  the  shell,  at  azimuth 
angles  of  (f,  90",  and  180"  (which  correspond  to  the  points  (0,0,1),  (0,1,0),  and 
(0,0,- 1),  using  the  coordinate  system  shown  in  Figure  1)  are  shown  in  Figure  8. 
Huang’s  new  70  term  Cesaro  sum  solution  [Ref.  14]  for  these  same  points  is  shown 
for  comparison  purposes.  While  our  finite  element  solution  shows  some  overshoot 
and  resulting  oscillation  for  0",  in  general  the  agreement  with  Huang’s  analytical 
solution  is  quite  good.  All  velocities  in  Figure  8  are  non-dimensionalized  to  be 
indep)endent  of  the  magnitude  of  the  incident  pressure  wave,  by  dividing  the 
original  non-dimensional  velocity  by  the  non-dimensional  magnitude  of  the  incident 
pressure  wave. 

B.  INFINITE  CYLINDER/PLANE  STEP  WAVE 

Figure  9  shows  the  geometry  of  the  infinite  cylinder/plane  step  wave 
problem.  The  same  material  properties,  parameters,  and  non-dimensionalization 
procedures  used  in  the  spherical  shell/plane  step  wave  problem  were  used  for  this 
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Figure  9.  Infinite  Cylinder/Plane  Step  Wave  Problem  Geometry 

problem,  except  that  a  shell  thickness  to  radius  ratio  of  0.029056  was  used  for  this 
problem. 

Because  of  the  symmetry  of  the  problem,  only  a  single  0.1  cylinder  radii 
wide  "ring"  of  the  infinite  cylinder  was  modeled  for  our  finite  element  analysis. 
In  addition,  since  the  problem  has  symmetry  about  the  plane  defined  by  a  point  on 
the  axis  of  the  cylinder  and  the  vector  normal  to  the  incoming  pressure  wave  front, 
only  one  half  of  this  ring  was  modeled. 

The  shell  was  modeled  with  36  Lagrangian  elements,  with  appropriate 
translational  and  rotational  symmetry  constraints  placed  upon  the  grid  points 
associated  with  these  elements.  MSC/DYTRAN’s  'General  Coupling  fluid- 
structure  interaction  method  was  used  for  this  problem;  to  form  the  closed  volume 
coupling  surface  required  for  this  method,  72  dummy  triangular  and  two  dummy 
quadrilateral  elements  were  defined.  Figure  10  shows  the  Lagrangian  structural 
elements  and  dummy  elements  used  in  our  finite  element  model  for  this  problem. 

The  fluid  mesh  used  for  this  problem  consisted  of  a  thin  block  of  elements. 
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Figure  10.  Structural  and  Dummy  Elements  in  Finite  Element  Model  for  Infinite 
Cylinder/Plane  Step  Wave  Problem 


with  the  dimensions  0.1  x  4  x  8  cylinder  radii.  This  fluid  block  was  meshed  with 
1  X  88  X  176  hexahedron  elements,  for  a  total  of  15488  fluid  elements.  Thus,  the 
length  of  each  element  in  the  y  and  z  directions  was  1/1 1  cylinder  radii.  The  fluid 
mesh  used  is  shown  in  Figure  11.  As  in  the  Spherical  Shell/Plane  Step  Wave 
problem,  the  amount  of  fluid  modeled  is  sufficient  to  prevent  reflection  of  acoustic 
waves  from  the  boundaries  from  affecting  the  solution  for  times  less  than  6  radius 
transit  times. 

All  boundary  conditions  for  the  fluid  mesh  shown  in  Figure  1 1  were  left  as 
default  "wall"  (no  flow)  boundaries,  except  for  the  boundary  at  z=4  cylinder  radii, 
which  was  given  a  "flow"  boundary  condition  with  a  prescribed  pressure  of  0.001 
bulk  modulii  and  a  z  direction  particle  velocity  of  -0.001  cylinder  radii/cylinder 
radius  transit  time.  Initial  conditions  were  prescribed  such  that  all  fluid  between 
the  z=4  and  z=l  cylinder  radii  planes  had  these  same  values,  and  the  remaining 
fluid  had  zero  initial  pressure  and  particle  velocity.  Thus,  time  t=0  corresponds  to 
the  instant  when  the  pressure  wave  just  touches  the  cylinder,  at  the  point  (0,0,1) 
(using  the  rectangular  coordinates  of  Figure  9;  in  cylindrical  coordinates  this  point 
is  at  a  radius  of  one  cylinder  radii  from  the  cylinder  axis,  at  an  angle  of  0^). 

The  size  and  position  of  the  cylinder  relative  to  the  fluid  volume  modeled 
is  illustrated  in  Figure  12.  Only  the  outline  of  the  fluid  volume  modeled  is  shown 
in  this  figure. 

The  entire  finite  element  model  for  this  problem  was  created  using  MSC/XL. 
The  MSC/XL  input  data  stream  and  MSC/DYTRAN  run  file  used  for  this  analysis 
are  provided  in  Appendix  B.  These  files  are  sufficient  to  duplicate  our  analysis. 

Results  from  our  analysis,  for  the  pressure-independent  non-dimensional 
radial  velocity  of  the  shell  at  0",  90",  and  180"  are  compared  with  the  8-term  finite 
series  analytical  solution  found  by  Huang  [Ref.  11]  in  Figure  13;  again,  very  good 
agreement  between  the  finite  element  and  analytical  solutions  is  seen. 
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Figure  13.  Non-Dimensional  Incident  Pressure  Independent  Radial  Velocity 
Versus  Non-Dimensional  Time  for  Infinite  Cylinder/Plane  Step  Wave  Problem 
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IV.  MULTIMATERIAL  EULERIAN  ANALYSIS  OF  BUBBLES 


As  mentioned  in  the  introduction,  the  whipping  induced  on  nearby  marine 
and  submarine  structures  due  the  oscillation  of  the  bubble  produced  by  an 
underwater  explosion  is  of  considerable  interest,  as  the  damage  due  to  this 
phenomena  can  exceed  that  due  to  the  primary  shock  wave.  In  order  to  investigate 
the  feasibility  of  modeling  the  physics  of  underwater  explosion  gas  bubbles  using 
a  multimaterial  Eulerian  hydrocode  formulation,  several  finite  element  analyses 
were  conducted. 

A.  DEEP  SPHERICAL  BUBBLE 

The  first  bubble  related  finite  element  analysis  conducted  investigated  the 
detonation  of  a  deeply  submerged  explosive  charge  in  the  absence  of  nearby 
boundary  surfaces.  As  deep  charges  are  known  to  undergo  little  vertical  migration 
due  to  gravity  [Ref.  15],  a  quasi  one-dimensional  finite  element  model,  in  which 
gravity  is  neglected  entirely,  was  used  for  this  analysis. 

The  problem  thus  has  spherical  symmetry.  The  particular  parameters  chosen 
for  this  analysis  correspond  to  one  of  a  series  of  tests  conducted  during  and  shortly 
after  World  War  II  at  the  Woods  Hole  Oceanographic  Institute  [Ref.  16],  and 
consist  of  a  0.299  kg  (0.660  lb)  TNT  charge  detonated  at  a  depth  of  178.6  m  (586 
ft).  The  geometry  of  this  problem  is  shown  in  Figure  14. 

The  TNT  in  this  problem  was  modeled  using  a  Jones-Wilkins-Lee  (JWL) 
equation  of  state,  with  standard  equation  of  state  and  detonation  velocity  parameters 
for  TNT  [Ref.  17].  With  this  state  equation,  the  pressure  in  the  "burned  fraction" 
of  the  explosive  material  is  related  to  the  specific  internal  energy  and  the  density 
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Figure  14.  Deep  Spherical  Bubble  Problem  Geometry 


by; 


e 


’I  +B|1 


_  WT]  '' 


R. 


1  ) 


h 

’I  +o)tiPo£’ 


(2) 


where 

n  =  p/po 

Po  =  initial  density 

E  =  specific  internal  energy  (per  unit  mass) 

A,B,a),Ri,R2  are  constants  for  the  explosive 

The  "burned  fraction"  is  just  that  portion  of  the  explosive  contained  within  the 
sphere  formed  by  the  spherical  detonation  wave  front  propagating  outward  from 
the  detonation  point  with  detonation  velocity  d.  As  the  nominal  0.660  lb  TNT 
charge  size  reported  for  this  experiment  represents  a  slight  increase  over  the  actual 
charge  size,  as  calculated  by  the  original  researchers  to  account  for  the  increased 
energy  of  the  booster  and  detonator  [Ref.  16],  the  explosive  was  modeled  as  a 
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homogeneous  mass  of  TNT.  The  equation  of  state  parameters  used  were  [Ref.  17]: 

Po  =  1630  kg/m"* 

E  =  4.29  X  10"  J/kg 

A  =  3.712  X  10“  Pa 

B  =  3.231  X  10^  Pa 

a>  =  0.30 

R,  =  4.15 

R,  =  0.95 

d  =  6390  m/s 

In  order  to  model  the  seawater  in  which  this  experiment  was  conducted,  a 
polynomial  equation  of  state  was  used.  This  state  equation  relates  the  pressure  in 
the  fluid  to  the  acoustic  condensation  g  and  the  specific  internal  energy  by: 


p=aiP+a2|J.^+a3P^+ Po-^  (M>0) 

(3) 

p=a^\i+  {bQ+b^\i)  PqE  (g<0) 

where 

p  =  (p  -  PoVpo 
Po  =  initial  density 

E  =  specific  internal  energy  (per  unit  mass) 
a,,a2,a3,bo,b|,b2  are  constants  for  the  fluid 

and  the  upper  equation  applies  to  a  fluid  in  a  compressed  state,  while  the  lower 

applies  to  a  fluid  in  an  expanded  state.  Constants  for  this  state  equation  were 

determined  by  using  the  mathematical  calculation  program  MATHCAD  [Ref.  18] 

to  fit  available  literature  equation  of  slate  data  for  fresh  water  compressed  to 

extreme  pressures  to  the  form  of  equation  (3),  and  then  replacing  the  initial 

(reference)  density  and  the  first  (adiabatic  bulk  modulus)  term  with  the  reference 

density  and  adiabatic  bulk  modulus  of  seawater.  This  determination  of  polynomial 

equation  of  state  constants  for  seawater  is  provided  in  Appendix  C.  The  seawater 

equation  of  state  parameter  values  determined  by  this  procedure,  appropriate  for 
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condensation  values  on  the  order  of  p<0.8,  are: 

ai  =  2.306  X  10®  Pa 

a,  =  8.432  X  10®  Pa 

aj  =  8.014  X  10®  Pa 

bo  =  0.4934 

bj  =  1.3937 

b2  =  0.0000 

Po  =  1025  kg/m' 

The  initial  conditions  given  to  the  seawater  modeled  in  this  problem  were 
an  initial  acoustic  condensation  of  zero  (initial  density  of  1025  kg/m')  and  an  initial 
specific  internal  energy  of  3750.4  J/kg.  This  initial  specific  internal  energy  value 
was  determined  from  equation  (3)and  the  above  constants;  it  represents  the  specific 
internal  energy  necessary  to  give  the  seawater  an  initial  pressure  equal  to  standard 
atmospheric  pressure  plus  the  pgh  pressure  at  a  depth  of  178.6  m  in  seawater,  using 
the  standard  gravitational  constant  value  for  g  of  9.80665  m/s . 

As  MSC/DYTRAN  is  a  three-dimensional  code,  the  increasing  volume  as 
one  moves  away  from  the  center  of  the  charge  was  taken  into  account  in  our  quasi 
one-dimensional  model  by  modeling  a  diverging  pyramid  shaped  region  of  fluid, 
with  default  "wall"  (no  flow)  boundary  conditions  on  all  sides.  However,  the  very 
end  of  this  pyramid  shaped  volume  was  modeled  using  rectangular  parallelepiped 
shaped  hexahedron  elements  with  the  same  volume  as  the  equivalent  pyramid. 
This  both  avoids  the  use  of  a  five  sided  pyramid  shaped  volume  at  the  point  of  the 
pyramid  (which  is  not  a  standard  finite  element  shape)  and  keeps  the  stable  time 
step  size,  which  by  the  Courant  Criterion  is  proportional  to  the  smallest  dimension 
of  an  element,  from  being  very  small.  Rectangular  parallelepiped  shaped 
hexahedron  elements  were  used  to  model  the  charge,  and  non-rectangular 
hexahedron  elements,  whose  top  and  bottom  were  parallel  but  whose  sides  had  a 
slope  of  0.1,  were  used  to  model  the  fluid.  As  our  primary  interest  was  in  the 
bubble  rather  than  the  initial  shock  wave,  only  three  elements  were  used  to  model 
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the  charge.  In  order  to  prevent  reflection  from  the  far  end  of  the  finite  volume  of 
fluid  modeled  for  the  duration  of  the  analysis,  a  large  volume  of  fluid  was 
modeled.  The  radial  dimension  of  the  seawater  in  the  model  extended  to  10,000 
times  the  initial  3.526  cm  "radius"  (height  of  rectangular  parallelepiped  actually 
modeled  was  1/3  this  radius)  of  the  charge.  To  model  this  large  volume  of 
seawater,  996  hexahedron  elements  with  a  non-uniform  radial  mesh  spacing  were 
used,  with  elements  becoming  larger  further  away  from  the  charge  center  (where 
gradients  are  smaller).  The  finite  element  model  used  in  this  analysis  is  shown 
from  a  three-dimensional  perspective  in  Figure  15. 

To  determine  the  radius  versus  time  behavior  of  the  explosion  gas  bubble, 
the  volume  of  explosive  gases  at  each  time  step  was  output.  This  volume  was  then 
converted  to  a  full  spherical  volume  by  multiplying  by  the  appropriate  geometric 
factor  and  then  used  to  calculate  the  radius  of  this  sphere.  The  resulting  radius 
versus  time  curve  is  shown  in  Figure  16.  Also  shown  in  this  figure  are  the 
experimentally  determined  first  maximum  bubble  radius  and  first  bubble  oscillation 
period  [Ref.  16].  The  finite  element  analysis  and  experimental  results  for  these 
quantities  are  seen  to  be  in  excellent  agreement. 

As  only  maximum  bubble  radii  and  bubble  periods  were  measured  in  this 
set  of  experiments,  it  is  possible  that  the  finite  element  analysis  radius  versus  time 
curve  could  agree  with  experimental  results  at  these  two  points  but  still  have  the 
wrong  shape.  However,  while  an  analytical  solution  does  not  exist  for  this 
problem,  we  can  calculate  an  analytical  solution  for  the  simplified  problem  of  a 
bubble  with  negligible  internal  energy  expanding  and  collapsing  in  an 
incompressible  media.  Expressing  the  total  energy  associated  with  radial  flow  of 
the  fluid  in  terms  of  the  maximum  radius  of  the  bubble  (mathematically,  this  is 
equivalent  to  using  the  maximum  radius  of  the  bubble  to  determine  the  initial 
conditions  of  the  problem)  and  then  separating  variables  gives  the  relationship 
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Figure  16.  Radius  Versus  Time  Behavior  for  Deep  Spherical  Bubble  Problem 
(0.2990  kg  TNT  Charge  at  178.6  m  Depth  in  Seawater) 


between  the  bubble  radius  and  time  as  [Ref.  19]: 


The  radius  versus  time  curve  obtained  by  numerically  integrating  this  equation  step- 
by-step,  using  the  seawater  density  and  hydrostatic  pressure  used  in  the  finite 
element  analysis  for  this  problem,  the  maximum  radius  obtained  from  that  analysis, 
and  an  initial  radius  of  zero,  is  also  plotted  in  Figure  16. 

Figure  16  thus  shows  that  not  only  are  the  maximum  bubble  radius  and 
bubble  period  obtained  from  the  finite  element  analysis  for  the  first  bubble 
oscillation  in  excellent  agreement  with  the  experimental  values  for  these  quantities, 
but  the  shape  of  the  radius  versus  time  curve  is  also  very  similar  to  that  obtained 
from  analytical  analysis  of  the  simplified  (incompressible/negligible  internal  energy) 
problem.  Up  to  the  time  of  the  first  maximum  radius,  these  curves  nearly  overlay 
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one  another,  while  after  that  time  the  finite  element  solution  begins  lagging  the 
simplified  problem  analytical  solution. 

The  effects  of  a  more  realistic  state  equation  for  the  seawater  and  the 
explosion  bubble  gases  (as  opposed  to  the  seawater  being  incompressible  and  the 
internal  energy  of  the  bubble  gases  being  zero)  are  clearly  shown  in  Figure  16.  In 
the  more  realistic  (finite  element)  analysis,  the  bubble  oscillation  is  asymmetric 
about  the  maximum  radius  and  the  bubble  period  is  lengthened.  An  approximate 
analysis  by  Herring  [Ref.  19]  shows  that  a  "first-order"  effect  of  compressibility  is 
to  cause  the  compression  portion  of  the  oscillation  cycle  to  take  longer  than  the 
expansion  portion  of  the  cycle,  which  is  what  was  observed  in  the  finite  element 
analysis.  In  addition,  the  radius  versus  time  behavior  for  the  finite  element  analysis 
shows  a  decreasing  maximum  radius  and  period  and  increasing  minimum  radius  for 
subsequent  oscillations.  This  is  due  to  the  energy  loss  as.sociated  with  radiation  of 
an  outward  traveling  pressure  wave  (the  "bubble  pulse")  when  the  bubble 
approaches  a  radius  minimum  [Ref.  19]. 

Up  until  the  time  of  the  first  minimum,  the  results  from  the  finite  element 
analysis  are  in  excellent  agreement  with  both  experimental  and  approximate 
analytical  results.  As  there  is  at  least  a  few  percent  uncertainty  in  the  detonation 
energy  for  any  particular  explosion,  the  agreement  with  experimental  results  is  as 
good  as  can  reasonably  be  achieved.  The  maximum  bubble  radius  and  bubble 
period  for  the  second  oscillation  cycle  were  also  measured  in  this  exj)eriment,  the 
measured  value  for  the  second  maximum  radius  was  11.6  inches  (about  29.5  cm), 
and  the  measured  time  for  the  second  minimum  was  30.85  msec.  There  thus 
apj)ears  to  have  been  an  energy  loss  when  the  bubble  was  near  its  nunimum  radius 
in  excess  of  that  attributable  to  acoustic  radiation.  It  is  believed  that  this  energy 
loss  was  due  to  heat  transfer.  As  Hicks  has  pointed  out,  deep  (non-migrating) 
bubbles  are  unstable  near  their  first  minimum,  and  photographs  show  that  numerous 
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at  this  time  numerous  needle-like  water  jets  spray  into  the  gas  bubble,  significantly 
cooling  the  hot  bubble  gases  [Ref.  20] .  The  energy  loss  in  the  gas  bubble  due  to 
this  phenomena  has  not  been  accounted  for  in  our  finite  element  analysis. 

The  input  and  postprocessing  files  used  in  this  analysis,  as  well  as  the 
FORTRAN  program  used  to  compute  the  solution  to  the  simplified  (incompressible 
media/negligible  internal  energy)  problem  are  provided  in  Appendix  D. 

B.  BUBBLES  NEAR  RIGID  BOUNDARIES 

A  number  of  interesting  phenomena  are  associated  with  the  expansion  and 
collapse  of  gas  bubbles  generated  by  underwater  explosions  in  the  vicinity  of 
boundaries.  These  include  migration  of  the  bubble  towards  or  away  from  the 
boundary  and  changes  in  the  shape  and  oscillation  period  of  the  bubble.  Each  of 
these  phenomena  can  effect  the  damage  done  by  the  bubble  on  marine  or 
submarine  structures. 

One  of  the  limiting  cases  is  that  in  which  the  boundary  behaves  as  a 
perfectly  rigid  wall.  This  case  is  of  interest  because  the  effect  of  the  wall  on  the 
bubble  approximates  the  effect  on  the  bubble  of  either  a  hard  bottom  or  a  nearby 
stiff  marine  or  submarine  structure.  The  simplest  situation,  in  which  the  rigid 
boundary  can  be  treated  as  an  infinite  plane,  is  appropriate  when  the  explosion 
takes  place  near  a  hard  bottom  or  when  the  explosion  takes  place  near  a  stiff 
structure  with  a  radius  of  curvature  that  is  large  compared  with  the  dimensions  of 
the  resultant  bubble. 

The  results  discussed  in  this  section  are  for  this  situation.  The  charge 
modeled  is  a  10.24  kg  cylinder  of  TNT,  20  cm  in  diameter  and  20  cm  high, 
submerged  in  seawater  at  a  depth  of  1000  m.  The  first  analysis  conducted  was  for 
the  "free-field"  case  in  which  the  rigid  boundary  is  absent  (rigid  boundary  at 
infinite  distance  from  center  of  charge);  this  analysis  yielded  a  maximum  free-field 
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Charge  (10.24  kg  TNT) 
Free-Field  Maximum  Radius 

P  (atm)  =  101325  Pa 
p  =  1025  kg/m^ 
g  =  9.80665  m/s^  I 
h  =  1000  m 

P^  =  P  (atm)  +  pgh 


Rigid  Wall  (various  standoffs  shown) 


Figure  17.  Geometry  of  Bubble  Near  Rigid  Boundary  Problems 

bubble  radius  of  70.98  cm.  Subsequent  finite  element  analyses  were  then 
conducted  with  a  rigid  wall  placed  at  4.015,  2.008,  1.374,  and  1.000  times  this 

distance  from  the  center  of  the  charge. 

The  geometry  of  these  problems  is  shown  in  Figure  17.  As  in  the  previous 
section  (the  "Deep  Spherical  Bubble"  problem),  gravitational  effects  on  the  motion 
of  the  bubble  were  neglected,  and  the  gravitational  constant  was  used  only  to 
determine  the  hydrostatic  pressure  in  the  seawater  surrounding  the  charge,  which 
is  assumed  to  have  a  constant  value  equal  to  the  hydrostatic  pressure  at  the  depth 
of  the  charge.  The  problem  is  then  axisymmetric,  with  the  axis  of  symmetry  being 
the  line  normal  to  the  rigid  wall  passing  through  the  axis  of  the  cylindrical  charge, 
the  problem  could  alternatively  be  made  axisymmetric  by  using  a  horizontal  rigid 
wall,  who.se  normal  is  the  gravitational  vector. 

The  diverging  geometry  of  the  problem  was  again  modeled  using 
geometrically  diverging  (wedge  shaped)  finite  element  models.  Quasi  two- 
dimensional  models  were  used.  Three  different  types  of  finite  element  models 
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were  used  in  the  analyses.  For  the  free-field  analysis,  a  quarter  symmetry 
axisymmetric  model  was  used;  half  symmetry  axisymmetric  models  were  used  for 
the  analyses  which  involved  a  nearby  rigid  wall.  The  half  symmetry  axisymmetric 
models  for  the  cases  in  which  the  rigid  wall  was  4.015, 2.008,  and  1.374  maximum 
free-field  radii  from  the  center  of  the  charge  were  modeled  using  only  Eulerian 
elements,  with  the  rigid  wall  being  a  simple  no-flow  boundary  condition;  the  fluid 
media  is  thus  modeled  only  on  one  side  of  the  wall.  The  case  in  which  the  rigid 
wall  is  at  1  .OCX)  free-field  radii  from  the  center  of  the  charge  was  modeled  by  using 
a  full  half  symmetric  axisymmetric  model  and  a  single  large  Lagrangian  steel  solid 
element,  with  interaction  between  the  Eulerian  and  Lagrangian  materials  being 
acomplished  with  MSC/DYTRAN’s  "General  Coupling"  method.  While  this  is 
somewhat  less  efficient  than  the  procedure  used  for  the  other  rigid  wall  cases,  it  is 
the  simplest  procedure  for  locating  the  rigid  wall  at  precisely  1 .000  free-field  radii 
from  the  center  of  the  charge. 

The  finite  element  mesh  for  these  analyses  consisted  of  three  different 
Eulerian  "regions."  A  refined  region  near  the  charge,  where  gradients  are  highest, 
was  used  to  ensure  reasonable  accuracy  in  capturing  the  motion  of  the  bubble. 
Outside  of  this  "central"  region  a  large,  less  refined  region  was  used  to  ensure  that 
reflection  from  the  non-rigid  wall  boundaries  did  not  effect  the  results  during  the 
duration  of  the  analyses.  This  "outside"  region  can  be  thought  of  as  "storing" 
energy  and  momentum  during  bubble  expansion,  which  then  effects  bubble  collapse 
after  the  momentum  field  reverses  direction.  A  "transition"  region  was  used  to 
connect  the  central  and  outside  regions. 

The  overall  geometry  of  the  finite  element  model  for  the  quarter  symmetry 
axisymmetric  free-field  analysis  is  shown  in  Figure  18,  from  a  three-dimensional 
perspective.  Figure  19  shows  the  complete  finite  element  model  for  this  analysis 
from  a  two-dimensional  perspective.  A  close  up  view  of  the  lower  left  comer  of 
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Figure  19.  Quarter  Symmetry  Axisymmetric  Finite  Element  Model  for  Free-Field 
Analysis 
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Figure  19,  which  includes  the  "central"  and  "transisition"  regions,  and  a  small 
portion  of  the  "outside"  region,  is  shown  in  Figure  20.  The  4  x  4  element  square 
block  (from  a  two-dimensional  perspective)  of  elements  in  the  lower  left  comer  of 
this  figure  are  initially  filled  with  TNT;  the  remaining  elements  are  initially  filled 
with  seawater.  The  axisymmetric  symmetry  axis  in  Figures  19  and  20  lies  along 
the  left  edge  in  these  figures.  The  quarter  symmetry  plane,  which  is  normal  to  and 
bisects  the  axis  of  the  charge,  lies  along  the  bottom  edge  in  these  figures.  The 
central  region  consists  of  9801  Eulerian  elements,  and  the  outer  region  contains 
13600  Eulerian  elements.  The  two-dimensional  height  and  width  of  the  elements 
in  the  central  region  is  0.025  m;  the  overall  central  region  dimensions  are  2.475  x 
2.475  m.  The  total  "wedge"  of  fluid  modeled  is  T  of  arc  of  a  solid  sphere  with 
a  radius  of  4(K)  m.  The  inverted-L  shaped  transition  region,  which  provides  a  3:1 
mesh  coarsening,  is  composed  of  269  Eulerian  elements. 

The  above  described  finite  element  model  was  used  as  the  basic  building 
block  for  the  models  used  in  the  other  analyses.  For  the  analysis  in  which  the  ngid 
wall  is  located  at  1  .(KK)  maximum  free  field  radius  from  the  center  of  the  charge, 
the  fluid  portion  of  the  model  was  created  by  reflecting  the  free-field  finite  element 
model  about  its  quarter  symmetry  plane.  The  rigid  boundary  was  then  created  by 
modeling  a  single  large  solid  steel  Lagrangian  element  and  using  MSC/DYTRAN  s 
General  Coupling  algorithmn.  The  finite  element  model  geometry  for  this  analysis 
is  shown  in  Figure  21,  from  a  three-dimensional  perspective.  Figure  22  shows  a 
close  up  view  of  the  -area  in  which  the  charge  is  located,  from  a  two-dimensional 
perspective;  the  shaded  elements  in  this  figure  initially  contain  TNT,  while  the 

remaining  elements  initially  contain  seawater. 

For  the  analyses  in  which  the  rigid  wall  was  located  at  1.374,  2.008,  and 
4.015  maximum  free-field  radii  from  the  center  of  the  charge,  fluid  elements  were 
added  to  the  area  between  the  free-field  quarter  symmetry  plane  and  the  rigid  wall 


34 


X 

A 


Figure  20.  Close  Up  View  of  Finite  Element  Model  in  Area  Near  Charge  for 
Free-Field  Analysis 
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Figure  22.  Close  Up  View  of  Finite  Element  Model  in  Area  Near  Charge  for 
Rigid  Wall  at  1  .(XX)  Free-Field  Radii  From  Center  of  Charge 
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plane.  The  wall  itself  was  then  modeled  by  just  making  the  boundary  condition 
on  the  element  faces  boardering  the  desired  rigid  wall  location  a  no-flow  boundary 
condition.  Figure  23  shows  the  finite  element  model  geometry  for  a  typical  case 
(rigid  wall  at  4.015  free-field  radii  from  center  of  charge),  from  a  three-dimensional 
perspective.  The  complete  finite  element  model  for  this  case  is  shown  in  Figure 
24,  from  a  two-dimensional  perspective.  A  close  up  view  of  this  model  in  the  area 
near  the  charge  is  shown  in  Figure  25.  The  finite  element  models  for  the  cases  in 
which  the  rigid  wall  is  located  at  1.374  and  2.008  free-field  radii  from  the  center 
of  the  charge  are  similar. 

The  same  state  equations  for  TNT  and  seawater  used  in  the  "Deep  Spherical 
Bubble"  analysis  were  used  for  the,se  analyses.  Detonation  was  initiated  at  the 
center  of  the  TNT  charge  (midpoint  of  the  axisymmetric  axis  of  the  cylindrical 
charge).  The  initial  density  and  specific  internal  energy  of  seawater  at  a  depth  of 
KKK)  m  was  1025  kg/m'  and  20076  Pa  m3/kg,  respectively;  these  values  when  used 
with  the  previous  equation  of  state  coefficients  give  the  correct  hydrostatic  pressure 
at  a  depth  of  1000  m. 

For  simplicity  in  model  generation,  the  two  angled  faces  in  the  models 
which  appear  to  meet  at  the  y  axis  (see  Figures  18,  21,  and  23)  are  actually  offset 
a  very  small  distance  (+/-  5  pm);  they  thus  make  an  angle  of  +!-  T  not  from  the 
xy  plane  at  z=0,  but  from  the  xy  planes  at  z  =  -t-5  pm  and  z=  -5  pm,  respectively. 
This  small  offset  cannot  be  seen  on  any  reasonable  scale  (it  is  small  even  in 
comparison  with  the  other  dimensions  of  the  smallest  elements),  and  permits  the 
model  to  be  meshed  using  only  hexahedron  elements.  Since  we  are  using  a  three 
dimensional  code  to  model  an  axisymmetric  two  dimensional  problem,  the  stable 
timestep  size  safety  factor,  normally  taken  as  about  0.666,  was  increase  by  a  factor 
of  about  28.6  (l/tan(r)),  the  difference  between  z-direction  and  x  or  y-direction 
element  lengths  for  elements  bordering  the  axisymmetric  symmetry  axis.  This 
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Figure  23.  Finite  Element  Model  Geometry  for  Rigid  Wall  at  4.0015  Free-Field 
Radii  From  Center  of  Charge 


39 


Figure  24.  Finite  Element  Model  for  Rigid  Wall  at  4.015  Free-Field  Radii  From 
Center  of  Charge 
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Figure  25.  Close  up  View  of  Finite  Element  Model  in  Area  Near  Charge  for 
Rigid  Wall  at  4.015  Free-Field  Radii  From  Center  of  Charge 


41 


immediately  reduces  the  solution  time  by  a  factor  of  28.6,  but  has  the  potential  to 
introduce  instability  if  the  fluid  in  elements  far  from  the  axis  have  a  much  higher 
combined  wave  and  partical  velocity  than  that  in  any  elements  bordering  the  axis. 
Such  a  situation  is  unlikely,  and  no  solution  instability  was  seen.  Run  times  of 
about  two  hours  on  a  RISC  based  Unix  workstation  were  obtained. 

Analysis  results  for  the  "volume-equivalent  spherical  radius"  (radius  of  a 
spherical  bubble  with  the  same  volume  as  the  actual  bubble)  versus  time  behavior 
are  plotted  in  Figure  26.  The  volume  equivalent  spherical  radius  versus  time 
behavior  for  the  free-field  case  and  the  cases  where  the  rigid  wall  is  located  at 
4.015  and  2.008  free-field  radii  (h*  =  4.015  and  h*  =  2.008)  from  the  center  of  the 
charge  is  shown  in  the  upper  graph,  while  the  lower  graph  shows  this  same 
information  for  the  cases  in  which  the  rigid  wall  is  located  at  1 .374  and  1 .000  free- 
field  radii  from  the  center  of  the  charge.  The  h*  =  2.(K)8  radius  versus  time  curve 
is  repeated  in  the  lower  graph  in  Figure  26  as  a  reference.  From  these  graphs,  the 
main  effects  of  the  rigid  wall  are  seen  to  be  to  increase  the  period  of  oscillation, 
to  reduce  the  maximum  volume,  and  to  increase  the  minimum  volume  and  the 
amount  of  time  when  the  bubble  is  near  minimum  volume  (to  spread  out  the 
minimum). 

An  approximate  analysis  by  Herring  [Ref.  19]  for  a  spherical  bubble  with 
negligible  internal  energy  expanding  and  collapsing  in  an  incompressible  fluid  in 
the  vicinity  of  a  rigid  wall  gives  the  modified  oscillation  period  as 


t'=j\ 


i+hzi. 

4h 


(5) 


where  T’  is  the  period  when  the  wall  is  present,  T  is  the  free-field  period,  h  is  the 
initial  distance  from  the  bubble  center  to  the  wall,  and  r^yg  is  the  average  radius  of 
the  bubble  over  one  oscillation  in  the  absence  of  the  wall.  Herring’s  analysis  thus 
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Figure  26.  Equivalent  Spherical  Radius  vs  Time  for  Explosion  Gas  Bubbles  Near 
a  Rigid  Wall 
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First  Period  (T*) 


Figure  27.  Non-Dimensional  First  Oscillation  Period  vs  (Non-Dimensional 
Standoff  Distance  From  Rigid  Wall)  ' 

predicts  that  the  bubble  period  will  increase  linearly  with  the  quantity  (1/h).  While 
our  analyses  predict  an  increase  in  the  bubble  oscillation  p)eriod  when  a  rigid  wall 
is  present,  results  from  our  analyses  indicate  that  the  increase  is  not  linearly 
proportional  to  the  inverse  of  the  standoff  distance  between  the  center  of  the  charge 
and  the  wall.  Results  from  our  analyses  for  T*  versus  1/h*,  where  T*  =  T’/T  is 
the  non-dimensional  period  and  1/h*  is  the  inverse  of  the  non-dimensional  standoff 
distance  (maximum  free-field  bubble  radius  divided  by  the  initial  standoff  distance 
between  the  center  of  the  charge  and  the  wall)  are  plotted  in  Figure  27.  The 
equivalent  non-dimensional  quantities  predicted  by  Herring’s  formula  (equation 
(5)),  where  the  quantity  ra^g  was  calculated  by  numerically  integrating  the  radius 
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versus  time  curve  from  our  free-field  anulysis  and  dividing  that  quantity  by  the 
period  determined  from  that  analysis,  is  shown  for  comparison  purposes. 

Examining  Figure  27,  it  is  seen  that  at  standoff  distances  less  than  about  two 
free-field  radii  (1/h*  greater  than  about  0.5)  the  periods  determined  from  our 
analyses  increase  almost  linearly  with  l/li*,  and  the  T*  versus  1/h*  curve  has  about 
the  same  slope  as  that  predicted  by  Herring’s  analysis.  However,  for  standoff 
distances  greater  than  about  two  free-field  radii  our  analyses  indicate  that  the  wall 
has  a  smaller  effect  than  that  predicted  by  Herring’s  analysis.  This  is  due  to  the 
compressibility  of  the  fluid  media  in  our  analyses.  In  fact,  because  of  the  finite 
wave  sjjeed  in  a  compresible  media,  there  is  a  finite  distance  beyond  which  the 
wall  can  have  no  effect  on  the  bubble  within  a  given  time  frame.  The  peak  sound 
speed  seen  during  our  analyses  was  3733  m/s;  in  conjunction  with  the  14.761  msec 
bubble  period  from  our  free-field  analysis,  the  calculated  standff  distance  beyond 
which  the  wall  can  have  no  effect  is  at  most  38.8  free-field  radii.  This  point  is 
indicated  by  the  second  open  circle  from  the  left  in  Figure  27  (and  in  other 
subsequent  graphs).  The  remaining  open  circles  each  represent  a  data  point  from 
our  analyses. 

The  displacement  of  the  center  of  mass  of  the  bubble  as  a  function  of  time 
was  found  from  results  from  the  finite  element  analyses  by  numerically  integrating 
the  material  velocity  time  history  data.  Analysis  results  for  the  displacement  of  the 
center  of  mass  of  the  bubble  as  a  function  of  time  are  shown  in  Figure  28.  In  this 
figure  a  positive  displacement  is  a  displacement  away  from  the  rigid  wall,  while 
a  negative  displacement  is  a  displacement  towards  the  wall;  the  rigid  wall  can  thus 
be  visualized  as  lying  below  the  time  axis.  The  upper  graph  in  Figure  28  shows 
the  displacement  of  the  center  of  mass  of  the  bubble  for  the  cases  in  which  the 
rigid  wall  is  located  at  4.015  and  2.008  maximum  free-field  radii  from  the  initial 
position  of  the  center  of  the  charge,  and  the  lower  graph  shows  this  quantity  for 
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Figure  28.  Displacement  of  Center  of  Mass  vs  Time  for  Explosion  Gas  Bubbles 
Near  a  Rigid  Wall 


initial  standoff  distances  of  1 .374  and  1  .OCX)  maximum  free-field  radii.  The  general 
characteristic  of  these  curves  is  a  small  initial  displacement  away  from  the  wall 
during  the  earlier  portion  of  the  oscillation,  and  then  a  much  larger  motion  towards 
the  wall  as  the  bubble  nears  minimum  volume.  This  general  curve  shape  has  been 
explained  by  Cole  [Ref.  20]  in  terms  of  the  wall’s  interference  with  the  flow  during 
expansion  and  collapse  and  the  resultant  momentum  field  in  the  fluid: 

In  the  case  of  a  rigid  surface,  the  presence  of  the  boundary  interferes  with 
radial  flow  of  water,  whether  outward  or  inward,  near  a  spherical  surface  in  its 
vicinity.  Initially,  when  the  pressure  in  the  gas  is  in  excess  of  the  hydrostatic 
pressure,  the  water  on  the  bubble  surface  near  the  wall  is  less  readily  displaced,  and 
the  bubble  surface  moves  away  from  the  wall.  The  effect  is  relatively  small, 
however,  because  the  net  pressure  (in  excess  of  hydrostatic)  is  positive  for  a  short 
part  of  the  bubble  period  ,  and  the  bubble  is  small  during  this  time.  When  the 
pressure  falls  below  hydrostatic,  acceleration  of  the  flow  toward  the  bubble  surface 
does  not  occur  as  readily  on  the  side  toward  the  wall,  and  the  flow  must  be  such 
as  to  bring  the  surface  nearer  to  the  wall.  A  considerable  amount  of  momentum 
is  imparted  to  a  large  mass  of  water  this  way  when  the  bubble  is  large.  As  the 
bubble  contracts,  the  momentum  aquired  becomes  concentrated  in  a  smaller  mass 
of  water  near  the  bubble,  and  the  velocity  of  flow  in  this  region  increases.  The 
bubble  surface  must  then  move  toward  the  wall  with  increasing  speed  as  if  attracted 
to  it.  This  effect  is  so  much  larger  than  the  repulsion  when  the  pressure  exceeds 
hydrostatic  that  the  dominant  motion  is  an  apparent  attraction  increasing  the  bubble 
velocity  toward  the  wall  as  it  contracts,  even  though  the  momentum  of  the  flow  is 
decrea.sing  in  the  most  contracted  stages. 

An  approximate  analysis  by  Herring  [Ref.  21],  in  which  terms  of  order 
higher  than  l/h*'  are  neglected,  the  bubble  is  assumed  to  remain  spherical  and  to 
have  negligible  internal  energy,  the  fluid  is  treated  as  incompressible,  and  the 
migration  of  the  bubble  is  treated  as  a  small  correction  to  the  motion  of  the  bubble, 
gives  the  rate  of  migration  of  the  center  of  the  bubble  as 


^  ^  3  f  ^j-4/  c?T 

dt  4/2 2  dt  2h^T^Jo  \  dx  I 


(6) 
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where  r(t)  is  the  radius  at  time  t  in  the  absence  of  the  wall.  In  conjunction  with 
equation  (5),  this  equation  can  be  seperated  an  integrated  to  give  the  distance  to  the 
wall  as  a  function  of  time.  However,  results  obtained  from  doing  so  are  known  to 
be  in  poor  agreement  with  experimental  results,  as  this  equation  is  very  sensitive 
to  the  bubble  radius  versus  time  behavior  when  the  bubble  is  near  minimum  radius, 
which  is  where  the  assumptions  used  in  its  derivation  (spherical  bubble  with 
negligible  internal  energy,  incompressible  media)  are  least  accurate  [Ref.  20]. 
Nevertheless,  this  equation  does  contain  a  periodic  term  which  changes  sign  at  the 
maximum  bubble  radius,  and  a  monotonic  term  which  becomes  large  near  the 
minimum  radius,  which  again  is  in  general  agreement  with  the  finite  element 
results  shown  in  Figure  28.  These  curves  are  also  in  qualitative  agreement  with 
results  obtained  by  Campbell  [Ref.  22]  from  the  shallow  water  detonation  of 
blasting  caps  near  a  vertical  wall. 

The  curves  in  Figure  28  also  show  a  relatively  small  magnitude  higher 
frequency  oscillation  superimposed  upon  the  lower  frequency  migration  curves. 
We  had  not  seen  this  phenomena  discussed  in  the  literature,  and  had  not  anticipated 
its  occurance.  The  primary  clue  to  its  cause  is  the  fact  that  there  appears  to  be  a 
linear  relationship  between  the  time  the  high  frequency  oscillations  begin  and  the 
standoff  distance  to  the  wall.  In  fact,  it  happens  that  the  time  the  oscillations  start, 
in  milliseconds,  is  roughly  equal  to  the  standoff  distance,  in  free-field  radii,  e.g.  the 
oscillations  start  at  about  two  milliseconds  for  the  case  where  the  charge  center  is 
initially  two  free-field  radii  from  the  wall.  Using  the  maximum  free-field  radius 
of  .7098  m  and  an  acoustic  wave  velocity  of  1500  m/s,  the  calculated  time  for  an 
acoustic  wave  to  travel  from  the  center  of  the  charge  to  the  wall  and  back  for  this 
case  is  1 .89  milliseconds,  which  is  about  when  the  oscillations  start.  The  higher 
frequency  component  is  thus  apparently  due  to  impingement  of  the  reflection  of  the 
initial  shock  wave  from  the  boundary  on  the  bubble. 
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The  presence  of  this  higher  frequency  component  in  the  migration  curves 
obtained  from  our  finite  element  analysis  therefore  appears  to  be  just  the  realization 
of  a  real,  physical  phenomena.  The  fact  that  the  oscillation  period  of  this  high 
frequency  component  increases  as  the  standoff  distance  increases  can  then  be 
explained  by  noting  that  for  larger  standoff  distances  the  time  required  for  the 
reflected  primary  shock  wave  to  reach  the  bubble  is  longer,  hence  the  bubble  has 
expanded  to  a  larger  radius,  and  so  the  "envelopment  time"  ( the  time  required  for 
this  shock  wave  to  transit  the  diameter  of  the  bubble)  is  longer. 

Other  characteristics  of  the  bubble  can  be  examined  in  terms  of  the  standoff 
distance  to  the  rigid  wall.  Figure  29  shows  the  maximum  and  minimum  spherical 
equivalent  bubble  radii  as  a  function  of  the  inverse  of  the  standoff  distance  (in  free- 
field  maximum  radii).  The  term  "equivalent  is  again  used  to  indicate  that  these 
are  the  radii  of  spherical  bubbles  with  the  same  volume  as  the  actual  bubbles,  the 
plotted  maximum  and  minimum  equivalent  radii  are  in  terms  of  the  maximum  and 
minimum  radii  for  the  free-field  bubble  (70.98  cm  and  22.71  cm,  respectively). 
Note  that  the  graphs  in  this  figure  have  different  vertical  scales,  so  that  while  the 
largest  reduction  in  the  maximum  equivalent  radius  is  less  than  5%  of  the  free-field 
value,  the  largest  increase  in  the  minimum  equivalent  radius  is  over  20%  of  the 
free-field  value.  By  analogy  with  the  free-field  case,  these  quantities  can  be 
viewed  as  being  related  to  the  energy  loss  of  the  system,  a  bubble  with  a  smaller 
maximum  radius  or  a  larger  minimum  radius  having  less  energy  available  for 
subsequent  oscillations. 

The  displacement  of  the  bubbles  at  the  end  of  the  first  oscillation  is  another 
quantity  of  interest.  The  bubble  pulse  emitted  when  a  bubble  is  near  its  minimum 
radius  can  be  approximated  as  eminating  from  the  location  of  the  center  of  the 
bubble  at  this  time,  and  if  appreciable  bubble  migration  has  occured  the  geometry 
of  the  bubble  pulse  loading  problem  is  then  changed.  The  top  graph  in  Figure  30 
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Displacement  of  C.M.  at  First  Minimum  Radius  (cm) 


shows,  for  the  cases  analyzed,  the  displacement  of  the  center  of  mass  of  the 
bubbles  at  the  the  time  of  their  first  minimum  radius,  plotted  as  a  function  of  the 
inverse  standoff  distance  (in  free-field  maximum  radii).  The  lower  graph  in  this 
figure  shows  the  peak  velocity  of  the  bubble  center  of  mass,  again  plotted  against 

the  non-dimensional  inverse  standoff  distance. 

The  shape  of  the  bubbles  for  the  various  cases  analyzed  was  also  examined. 
In  order  to  provide  for  meaningful  comparison  between  the  cases,  the  bubble  shape 
was  plotted  at  the  same  fractions  of  the  first  oscillation  periods  rather  than  at  the 
same  actual  times.  As  the  first  oscillation  period  is  different  for  each  standoff 
distance,  a  second  finite  element  analysis  of  each  model  was  done.  The  first  period 
for  each  case  was  found  from  the  first  analysis,  which  was  then  used  to  determine 
the  times  to  ask  for  a  complete  dump  of  the  element  densities  during  the  second 
analysis.  Plots  of  the  bubble  shapes  are  provided  in  Appendix  E,  at  times  of 
0.5,0.8,().9,0.95,().98,().99,  and  1.0  first  bubble  periods.  These  plots  show  that  for 
standoff  distance  of  about  two  free-field  radii  or  more,  the  bubble  remains  fairly 
spherical  up  to  the  time  of  the  first  minimum.  In  the  cases  for  which  the  rigid  wall 
is  initially  closer  than  two  free-field  radii,  as  the  bubble  approaches  minimum 
volume  it  assumes  a  "kidney"  shape,  similar  to  that  seen  experimentally  for  bubbles 
migrating  due  to  gravity  [Ref.  23].  For  the  case  in  which  the  rigid  wall  is  initially 
one  maximum  free-field  radii  from  the  center  of  the  charge,  the  kidney  shape 
evolved  into  a  torus,  indicating  penetration  of  a  wall-directed  jet  through  the  center 

of  the  bubble. 

The  input  and  postprocessing  files  used  for  the  cases  analyzed  in  this  section 
are  provided  in  Appendix  F. 
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V.  CONCLUSION 


This  report  describes  analysis  procedures  used  in  and  results  obtained  by 
directly  applying  coupled  Lagrangian-Eulerian  and  multimaterial  Eulerian  finite 
element  analysis  to  several  problems  of  interest  in  underwater  shock  research.  The 
problem  types  analyzed  encompass  both  classical  acoustic  wave-shell  fluid-structure 
interaction  and  analysis  of  bubbles  produced  by  the  detonation  of  explosives 
underwater. 

The  results  for  the  acoustic  wave-shell  fluid-structure  interaction  problems 
compared  quite  well  with  the  analytical  solutions  for  these  problems.  Results  for 
the  bubble  problems  analyzed  were  also  quite  satisfactory.  For  the  one-dimensional 
free-field  (spherical)  bubble,  results  for  the  first  bubble  maximum  radius  and  period 
were  in  very  good  agreement  with  results  from  the  experiment  being  modeled,  and 
the  shape  of  the  bubble  radius  versus  time  curve  was  similar  to  the  theoretical 
result  for  a  bubble  with  negligible  internal  energy  in  an  incompressible  media.  For 
the  other  bubble  problem  analyzed,  an  explosion  generated  gas  bubble  at  various 
distances  from  a  rigid  wall,  there  were  fewer  theoretical/experimental  results 
available  with  which  to  make  meaningful  comparisons;  the  results  obtained  were 
in  general  agreement  with  those  that  were  available.  The  large  amount  of  data 
available  from  the  analyses  of  a  bubble  at  various  initial  standoff  distances  from 
a  rigid  wall  provide  a  useful  characterization  of  the  effects  of  the  standoff  distance 
on  the  bubble. 

One  of  the  benefits  of  the  direct  finite  element  method  used  is  that  it  does 
not  rely  on  time  or  frequency  domain  approximations,  so  that  the  solution  accuracy 
obtained  is  dependent  only  upon  the  fineness  of  the  mesh  used  and  the  accuracy 
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with  which  the  state  equation  parameters  for  the  materials  modeled  are  known. 
This  method  does  involve  a  far  greater  number  of  elements  than  a  boundary 
element  method,  however,  the  ever  expanding  capability  of  computers  makes  direct 
application  of  the  finite  element  method  using  coupled  Lagrangian-Eulerian  and 
multimaterial  Eulerian  analysis  practical  for  an  increasing  number  of  problems. 
Reasonable  problem  solution  times  can  be  obtained,  since  time-marching  using  an 
explicit  finite  difference  technique  can  be  very  efficient  even  for  very  large 
problems  (no  eigenvalue  problem  need  be  solved). 
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APPENDIX  A:  SPHERICAL  SHELL/PLANE  STEP  WAVE  INPUT 


$  MSC/DYTRAN  RUK  FILE  :  Spherical  Shell/Plane  Step  Wave 
$ 

START 

TIME  =  9999 
LIMMEM  =  100000 
CEND 
$ 

TITLE  =  Spherical  Shell /Plane  Step  Wave 

$ 

ENDSTEP  =  99999 
ENDTIME  =  6. 

PARAM,INISTEPi l.E-2 
PARAM , MINSTEP , 1 . E “ 9 
$ 

TYPE  (NODES)  =  TIMEHIS 
GRIDS  (NODES)  =  1 
SET  1  =  1,37,97 
GPOUT  (NODES)  =  YVEL,ZVEL 
STEPS  (NODES)  =  0, THRU, END, BY, 1 
SAVE  (NODES)  =  99999 
$ 

TLOAD  =  1 
SPC  =  1 
$ 

BEGIN  BULK 

$ 

SETl, 99,1, THRU, 600 
SURFACE, 100,, ELEM, 99 
COUPLE, 1, 100, INSIDE, ON, ON 
$ 

TICEUL  1 

+  CYLINDERl  2 

+  CYLINDER2  2 

CYLINDERl  0 . 0 

+  10.0 

CYLINDER2  0 • 0 

+  10.0 

TICVAL, 1, , ZVEL, - . 001 , PRESSURE , . 001, DENSITY, 1 . 001 
TICVAL , 2 , , ZVEL ,0.0, PRESSURE ,0.0, DENSITY ,1.0 
$ 

INCLUDE  buDc.dat 

$ 

ENDDATA 

$ 


1 

2 

0.0 

0.0 


9. 

1. 

5.0 

5.0 


0.0 

0.0 


0.0 

0.0 


+ 

+ 

1.0  + 

-5.0  + 
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c 

c  LS-INGRID  INPUT  FILE  :  Spherical  Shell/Plane  Step  Wave 
c 


c  This  file  is  used  to  create  all  three  "parts'*  for  this  problem 
c  by  "commenting  out"  2  of  the  3  parts  with  {curly  braces}  so  that 
c  only  1  part  is  read  in;  this  is  repeated  during  execution  of 
c  LS“ INGRID  for  each  part  in  turn 
c 

dn3d  vec 

mat  1  type  3  ro  1  e  1  pr  . 3  sigy  1  shell  endmat 
mat  2  type  3  ro  1  e  1  pr  . 3  sigy  1  endmat 

c 

c  Part  1  :  Spherical  Shell 
c 
{ 

start 
-1  6  -11; 

-1  6  -11; 

-1  6  -11; 

-10  1 
-10  1 
-10  1 

sfi  -1  -3;  -1  -3;  -1  -3;  sp  0  0  0  1.0 

d  1  0  0  2  0  0 
d  0  1  0  0  2  0 
mate  1 
end 
} 

c 

c  Part  2  :  Duinmy  Elements 
c 

start 
-1  6  -11; 

-1  6  -11; 

-1  6  -11; 

-10  1 
-10  1 
-10  1 

8fi  -1  -3;  -1  -3;  -1  -3;  sp  0  0  0  1.0 

a  2  2  0  3  3  0 

mate  1 
end 

c 

c  Part  3  :  Fluid  Elements 
c 
{ 

start 

1  17  33; 

1  17  33; 

1  17  49  65; 

0  2  4 
0  2  4 
-4-224 
mate  2 
end 
} 

end 
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INGRID>NASTRAN  TRANSLATION  PROGRAM  :  Spherical  Shell/Rlane  Step  Wave 

this  program  converts  output  from  the  preprocessor 
a  MSC/NASTRAN  compatible  geometry  data  file  (only  one  material 
can  be  translated  in  this  version) 

program  ing2nast 

iitplicit  double  precision  (a“h,o-z) 
character's  a 

open(14/file='ingrido' , status = 'old' ) 

open ( 15, file=' geometry. dat status=' new' )  » 

write{*,*)  'input  desired  NASTRAN  property  number  (PID) . 

read(*,*)  npid  ,  , 

write {*,*)  'input  gridpoint  numbering  offset: 
read{*,*)  ngoffset 

writ0(*,*)  'input  element  numbering  offset:' 
read ( * , * )  neof  f  set 
read (14,*) 

read { 14 , * )  tempi , ngpts , nhexas , temp2 , nquaas 
do  100  i=l,26 
read(14,*) 
continue 

do  200  i=l, ngpts 

read(14,991)  npt,a,x,y,z 

npt  =  npt+ngoffset  ^ 

write(15, 998) 'GRID*  '.npt,'  ^O^T 

write (15, 999 ) '*CONT  ',z 
continue 

if  ( nquads . gt . 0 )  then 
do  300  i=l, nquads 

read (14, 993)  nel,nmat ,nl,n2,n3 ,n4 

read (14, *) 

nel=nel+neof f set 

nl  =  nl+ngoffset 

n2  =  n2+ngoff8et 

n3  =  n3+ngoffset 

n4  =  n4+ngoffBet  ^  ^  ^  ^  a 

write(15,994)  'CQUAD4  Snel,npid,nl,  n2  ,  n3 ,  n4 
cont inue 
else 

do  400  i=l, nhexas  ,  r  *7  q 

read (14,995)  nel , nmat , nl , n2 , n3 , n4 , n5 , n6 , nv , n8 
nel  =  nel+neoffset 
nl  =  nl+ngoffset 
n2  =  n2+ngoffset 
n3  =  n3+ngoffset 
n4  =  n4+ngoffset 
n5  =  n5+ngoffset 
n6  =  nS+ngoffset 
n7  =  n7+ngoffset 

n8  =  nS+ngoffset  xl 

write(15,996)  'CHEXA  ' ,nel,npid,nl,n2,n3,n4,n5,nb,  + 

write(15,997)  '+  ',n7,n8 

cont inue 
endif 

format (i8 , a5, 3e20 . 0) 
format (i8, i5,4i8) 
format (a8 , 6i8) 
format (i8, is, 8i8) 
format (a8, 8i8,al) 
format (a8 , 2i8 ) 

format (a8 , il6 , al6, 2el6 . 9 , a8) 

f ormat (a8 , el6 . 9 ) 

close ( 14) 

close (15) 

end 
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!  MSC/XL  INPUT  DATA  STREAM  :  Spherical  Shell/Plane  Step  Wave 
! 

Read  MSCInput  File="shellgeo.dat"  Of f 8etType=None  /NoExec/NoCase/Bulk 

I 

Select  Part/1 

Read  MSCInput  Fil©="dumnry’geo.dat "  Of f s©tType=None  /NoExec/NoCase/Bulk 

Refresh/Find  View/1  /Erase /NoCenter/Complete/WireFrame/Undeformed/NoPlot/NoArrow 

Check  Grid/ltl602  Type=Duplicates  View=l  PartList=Otl  TolerancesO . 005 
Equivalence  Grids  /Update /NoExtend/NoCollaps©  _ 

Define  Group  Name=xzsygrd0  Type^Grid  View=l  PartLi8t=0tl  WindowType=RubberBandWindow 
CollectMode=Inside  Boundary= Include  RubberbandXMin=0 . 050091779  RubberbandXMax=l . 0593952 
RubberbandYMin=~0 . 042142798  RubberbandYMax=0 . 026815818  _ 

Define  Group  Naine=yzsygrds  Type=Grid  View=l  PartList=0tl  WindowType=RubberBandWindow 
Co llectMode= Inside  Boundary = Include  RubberbandXMin=-0 . 02200132  RubberbandXMaxsO . 043822813 
RubberbandYMin=0 . 036219266  Rubb©rband‘yMax=l .  0455226 

Define  Group  Nain€=zaxi8sygrds  Type=Grid  Vi©w=l  PartLi8t=0tl  WindowType=RubberBandWindow 
Coll©ctMode= Inside  Boundary=Include  RubberbandXMin=:-0 . 037673733  RxibberbandXMax^O  .  059495226 
RubberbandYMin= -0.051546246  RubberbandYMax=0 . 029950301 

Edit  SPC/1  GridList=xz8ygrds  DOF="246"  /Create/Update/NoExtend/NoOverwrite 
Edit  SPC/1  GridListsyzsygrds  DOF="156’*  /Create/Update/NoExtend/NoOverwrite 
Edit  SPC/1  GridList=zaxis8ygrds  DOF="12456"  /Create /Update/NoExt end/ NoOverwr it e 

I 

Select  Part/2 

Read  MSCInput  File=" fluidgeo.dat"  Of f setType=None  /NoExec/NoCase/Bulk 

Refresh/Find  View/2  /Erase/NoCenter/Complete/WireFrame/Undeformed/NoPlot/NoArrow 

Define  Group  Naine=leftbdryel8  Type=Element  View=2  PartList=0t2  WindowType=RubberBandWindow 
CollectMode=Inside  Bounda2ry= Include  RubberbandXMin=-41788 . 899  RubberbandXMax=-39224 . 691 
RubberbandYMin=-0 . 11872614  RubberbandyMax=4 . 2404263 

Define  CFace  ElementList  =  lef tbdryele  FaceNo=:4  FEPaceSet=l  PartLi8t=0t2 
OutputLi9t=1537t2560 

Edit  Flow/1  LoadSet=l  FEFaceSet=l  ZVel=-0.001  Pressure=l . 001  Density=l . 001 
Edit  TLoadl/1  SID=1  LID=1  DynaType=Flow  TableId=Blank  /Create 

I 

Edit /Create  DMatep/1  DMatRHO=7 . 7885  DMatE=89.097  DMatNU=:0.3 
Edit/Create  PShell/1  MID1=1  T=0. 029056 

Edit /Create  PShelll/2  FORM=Dxiinmy  QUADPS=Gau8s  REF=Mid 

Edit /Create  EOSPol/1  EOSAl=l 

Edit /Create  DMat/2  DMatRHO=l  EOID^l 

Edit /Create  PEulerl/3  TYPE=Hydro  TICEulID=l 

I 

Write  MSCInput  Fil©="bulk.dat "  Fonnat="Bulk"  /NoExec/NoCase/Bulk 
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APPENDIX  B:  INFINITE  CYLINDER/PLANE  STEP  WAVE  INPUT 


$  MSC/DYTRAN  RUN  FILE  :  Infinite  Cylinder /Plane  Step  Wave 

$ 

START 

TIME  =  999 
CEND 
$ 

TITLE  =  ELASTIC  CYLINDER  /  STEP  WAVE 

$ 

ENDSTEP  =  9999 
ENDTIME  =  6. 

PARAM, INISTEP, 1. 64E-2 

$ 

TYPE  (NODES)  =  TIMEHIS 
GRIDS  (NODES)  =  1 
SET  1  =  1,37,73 
GPOUT  (NODES)  =  ZVEL,YVEL 
STEPS  (NODES)  =  0 , THRU, END, BY, 1 
SAVE  (NODES)  =  9999 
$ 

TLOAD  =  1 
SPC  =  1 
$ 

BEGIN  BULK 

$ 

SETl, 3,1, THRU, 110 
SURFACE,!, ,ELEM,3 
COUPLE, 1, 1, INS IDE, ON, ON 
$ 

INCLUDE  bulk.dat 

$ 

ENDDATA 

$ 
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I 

i  MSC/XL  INPUT  FILE  :  Infinite  Cylinder/Plane  Step  Wave 

I 

Define  Point  X=-0.05  Y-0  Z=1 

Translate  Point/1  DeltaX=0.1  DeltaY=0  DeltaZ=0  /Create 
Translate  Point/lt2  DeltaX=0  DeltaY=:0  D6ltaZ=-l  /Create 

Connect  Points  PointlList=:l  Point2Li8t=2  «  ,  i  c  ^  n 

Rotate  Curve/1  FroinX=0  FromY=0  FromZ=0  ToX=l  Toy=0  ToZ=0  Angle=-5  Of f setAngle-O  /Create 
Connect  Curves  CurvelList=l  Curve2List=2  /Points  „  «  ^  .  r  «-n 

Rotate  35  Surface/1  FroinX=0  FromY=0  FromZ=0  ToX=l  ToY=0  ToZ=0  Angle=-5  OffsetAngle- 

/Create  ^  ^ 

MeshParetm  Surface/lt36  Type=Quad4  U=1  V=1  Pattern=l  PID=1 
Select  Part/1 

Connect  Points  PointlList=l,  2  Point2Li8t=3 , 4  ^  ,  e  yr-r<aiahe 

Rotate  Curve/3t4  FromX=0  FromY=0  FroinZ=0  ToX=l  ToY=0  ToZ=0  Angle=-5  OffsetAngle-  / 
Connect  Curves  CurvelList=:3 , 4  Curve2List=5 , 6  «  .  ,  c 

Rotate  35  Surface/37t38  FromX=0  FroinY=0  FromZ=0  ToX=l  ToY=0  ToZ=0  Angle=-5  Of f setAngle_0 

/Create 

NeshParam  Surf ace/37tl08  Type=Tria3  U=1  V=1  Pattern=l  PID=2 
Connect  Curves  CurvelList=3  Curve2Li8t=4 

Translate  Surface/109  DeltaX=0  DeltaY=0  DeltaZ=-l  /Create 

MeshParam  Surf ace/109tll0  Type=Quad4  U=1  V=:l  Pattern=l  PID=2 

Check  Grid/lt368  Type=Duplicate8  View=l  PartLi8t=0  Tolerance=0 . 005  /Update 

Equivalence  Grids  /Update /NoExtend/NoCollapse 

Renxunber  Grid/ltl49  PartList=0  OutputList=150t225  /Update 

Renumber  Grid/150t225  PartList=0  OutputList=lt76  /Update 

Edit  SPC/1  GridLi8t=l,2,73,74  DOF="12456"  /Create 

Edit  SPC/1  GridList=3t72  DOF="156"  /Create 

Select  Part/2 

Define  Point  X=-0.05  Y=0  Z=4  CID=0  n  i- 

Translate  Point/7  DeltaX=0  DeltaY=4  DeltaZ=0  CID=0  /Create 
Connect  Points  PointlList=7  Point2List=8 

Translate  Curve/7  DeltaX=0.1  DeltaY=0  DeltaZ=0  CID=0  /Create 
Connect  Curves  CurvelList=7  Curve2List=8  /Points 
Translate  Surface/111  DeltaX=0  DeltaY^O  DeltaZ=-8  CID=0  /Create 
Connect  Surfaces  Surf acelList=lll  Surf ace2List=:112 

MeshParam  Solid/1  Type=Hexa  0=88  V=1  ”=1’®  .  .  on 

Define  CFace  ElementList=llltl98  FaceNo=3  FEFaceSet=l  PartList-Otl 

Define  CFace  ElementLi0t=198tl5598b88  FaceNo=5  FEFaceSet=2  PartLi0t=Otl 

Define  CFace  BlementList=15511tl5598  FaceNo=6  FEFaceSet=2  PartLiBt=0tl 

Edit  Flow/1  LoadSet=l  FEFaceSet=l  ZVel=-0.001  Pressure=0 . 001  Denaity=l . 001  /Create 

Edit  Flow/2  LoadSet=l  FEFaceSet=2  /Create 

Edit  TLoadl/1  SID=1  LID=1  DynaType^Flow  /Create 

Edit  TICVal/1  Name= "DENSITY"  Value=l  Pair=l  /Create 

Edit  TICVal/1  Naine="ZVEL"  Value=0  Pair=2  /Modify 

Edit  TICVal/2  Name =" DENSITY"  Value=1.001  Pair=l  /Create 

0^4  4-  TTCVal/2  Name="ZVBL"  Values “0.001  Pair=2  /Modify  ,  ,  - 

Edit  TICEul/1  Type=Elements  TICEulLiBt=5919tl5598  MID=2  TICValId=l  Level=l  DMaterial-1 

EdirTICEul/l  Type=Elemente  TICEulLiat=lllt5918  MID=2  TICValId=2  Level=2  DMaterial=2 

Edit^DMatep/1  DMatRH0=7 .7885  DMatB=89.097  DMatNU=0.3  /Create 
Edit  PShell/1  MID1=1  T=0. 029056  /Create 

Edit  Pshelll/2  FORMsDummy  QUADPS=Gauss  REF=Mid  /Create 

Edit  EOSPol/1  EOSAlsl  /Create 

Edit  DMat/2  DMatRHO=l  EOID=l  /Create 

Edit  PEulerl/3  TYPE=Hydro  TICEulID=l  /Create  nv 

Write  MSCInput  Files  "bulk .  dat"  Foarmats"Bulk"  /NoExec/NoCase/Bulk 
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APPENDIX  C:  A  POLYNOMIAL  STATE  EQUATION  FOR  SEAWATER 


A  straightforw'ard  process  for  fitting  Giiineisen  state  equation  data  to  the  polynomial  state  equation 
usually  found  in  most  advanced  finite  element  programs  is  discussed,  and  polynomial  state  equation 
parameters  for  pure  water  and  seawater  are  obtained  from  published  values  of  the  Gruneisen  state  equation 
parameters  for  pure  water. 

GRUNEISEN  EQUATION  OF  STATE 

The  Gmneisen  state  equation  relates  the  pressure  pg.  the  condensation  |i,  and  the  specific  internal 
energy  E  by: 

In  compression: 


In  tension: 

Pg(F-F,ro.p  oX,aj  ;=p  o 

In  these  equations,  pg  is  the  local  pressure,  which  is  a  function  of  the  condensation  p  (p=(p*pQ)/pQ), 
related  by  the  six  parameters  7q,  Po^^,  a,  S  S2,  and  S3,  and  the  specific  internal  energy^  per  umt  volume  E 
and  condensation  p,  related  by  the  two  parameters  and  'a' .  Further  explanation  ot  this  state  equation  is 
provided  by  Gurtman  et  al  [Ref  24]  and  by  Hailquist  and  Stillman  [Ref  25]. 

POLYNOMIAL  EQUATION  OF  STATE 

A  typical  FEM  program  polynomial  state  equation  is: 

In  compression 

p  p(p,E,a  i,a  2. a  3,b  i,b  2)  pp-ha  2-P  +3  3'F  rP-fl>  2'F  )-E 

In  tension 

Pp(p,E.a  pbQ.b  1)  :=a  i'P-p(bo-hb  i  p)-E 

In  these  equations,  pp  is  a  function  of  p.,  related  by  three  parameters  (a| ,  a2,  and  33),  and  E  and  p,  related 
by  the  parameters  bQ,  b| ,  and  b2- 
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CURVE  FITTING 


The  appropriate  coefficients  for  converting  Gruneisen  state  equation  data  to  polynomial  state  equation 
data  can  be  determined  as  follows: 

Matching  terms  multiplying  E  yields: 

bo  =  1(0 

bj  =  a 
b2  =  0 

Matching  the  tension  equations  yields: 

ai  =poC2 

The  remaining  polynomial  coeficients  33  and  a3  can  then  be  found  by  fitting  the  terms  not  multiplymg 
E  in  the  Gmneisen  state  equation  for  compression  with  the  teims  in  the  polynomial  state  equation  also 
not  multiplying  E,  over  the  range  for  which  the  Gruneisen  state  equation  is  valid  or  the  range  the 
polynomial  state  equation  is  to  be  used.  This  procedui'e  is  earned  out  below  tor  pure  water. 


WATER  EOS  FOR  0<p<0.8  (po<P<l-8po) 

Si:=2.56  S2--1.986  S;,:=0.2268 

a:  =  1.3937  yo:=0-4934 

Po:  =  1.0  c:=0,1415 


Gruneisen  Coefficients  (g  -  cm  -  ps) 
[Ref  26,  27] 


Then  letting  p  vaiy  discretely  over  the  range  of  validity  of  the  Gruneisen  state  equation, 
and  defining  the  function  g(p)  that  is  to  be  approximated  with  the  cubic  function  h(p)  (wfoc 
has  coefficients  a  [ ,  a2,  and  a3  -  the  terms  in  both  state  equations  multiplying  E  are  already 

matched),  we  can  plot  g(p)  as  a  function  of  p. 


n;=0.0.0.01..0.8 
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An  examination  of  this  plot  shows 
that  a  cubic  polynomial  should  do 
a  reasonable  job  of  approximating 
this  function 


0  0.2  0.4  0.6  0.8 

A  \’ei‘y  simple  method  to  determine  the  coefficients  32  and  a^  (a^  is  already  known)  such  that 
the  cubic  polynomial  h(p)  approximates  g(p)  in  this  case  is  to  match  the  functions  h  and  g  at 
two  non- zero  points.  The  points  p  =  0.3  and  p=0.6  are  a  reasonable  choice. 


gW 


a  1-0.02002 


p  y  :-0.3  p 


Solving  2-equations  in  2-unknowns  then  gives 


2 

2 

r2 


3 

Pi  1 

^2  I 


\g(^l  ^  2/ 


1^2]  /  0.08432  \ 
Wv  \  0.08014/ 


We  can  then  plot  h(U),  the  residual  error,  and  the  %  conversion  error: 
h(p)  :=a  i-p-pa2P^4-a3  P 


residual: 
%  error: 


r(p)  ‘^gW-hW 

sW:  =  100-44 
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5 


Examining  these  plots,  we  see  that  this  is  a  reasonably  good  conversion  of  the  Gruneisen 
state  equation  to  a  polynomial  equation  of  state,  particularly  for  typical  condensation  values 
of  less  than  0.5  .  It  could  obviously  be  "tuned"  for  less  error  over  a  smaller  condensation 
range  of  interest.  In  addition,  a  more  sophisticated  procedure  for  finding  the  best  values  for 
32  and  a^  could  be  used,  i.e.  minimize  the  square  of  the  residual  error  over  the  condensation 

range  of  interest  fonnally.  In  this  case,  the  small  difference  between  the  two  state  equations 
obtained  is  probably  already  within  the  margin  of  uncertainty  which  results  from  the 
uncertainty  with  which  the  parameters  in  the  Gruneisen  state  equation  are  known,  so  further 
refinement  was  not  carried  out. 


SEAWATER  EOS  FOR  0<p<0.8  (po<p<1.8po) 

To  extend  the  above  polynomial  state  equation  for  use  with  seawater,  the  density  and  acoustic 
wave  velocity  of  seawater  were  used  in  place  of  the  density  and  acoustic  wave  velocity  of  pure  watei . 
Taking  these  as  1.025  g/cm-’  and  0.1500  cm/ps,  the  a;  coefficient  becomes 


aj;=0.02306 

The  remaining  constants  determined  above  for  pure  water  were  used  for  seawater;  this  polynomial  EOS 
is  thus  a  first-order  correction  to  the  pure  water  polymomial  EOS. 
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APPENDIX  D:  DEEP  SPHERICAL  BUBBLE  INPUT 


$  MSC/DYTRAN  RUN  FILE  :  Deep  Spherical  Bubble 

$ 

START 
TIME  =  999 
CEND 
$ 

TITLE  =  SPHERICAL  BUBBLE 

$ 

ENDSTEP  =  99999 
ENDTIME  =0,05 
PARAM, INISTEP , 5 . E>7 
P ARAM , RHOCUT , 1 . E - 6 
PARAM , ROH YDRO , 1 . E - 6 
PARAM , ROMULTI , 1 . E - 6 
$ 

TYPE  (MATHIS)  =  TIMEHIS 
MATS  (MATHIS)  =  1 
SET  1=1 

MATOUT  (MATHIS)  =  VOLUME 
STEPS  (MATHIS)  =  0 , THRU, END, BY, 10 
SAVE  (MATHIS)  =  99999 
$ 

BEGIN  BULK 

$ 

INCLUDE  bulk.dat 

$ 

ENDDATA 

$ 
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!  MSC/XL  INPUT  FILE  :  Deep  Spherical  Bubble 

I 

Edit  ApplicationTable  DuplicateElementToleraiice=0 . 0001 
Edit  ApplicationTable  DuplicateGridTolerance=0 . 0001 
Edit  ApplicationTable  DuplicateFEFaceTolerance=0 . 0001 
Edit  ApplicationTable  Degenerates lement Tolerance =0 . 0001 

Define  Point  X=0. 03526  Y=-0. 001763  Z=-0, 001763  CID=0 

Translate  Point/1  DeltaX=0  DeltaY=0 . 003526  DeltaZ=0  CID=0  /Create 

Connect  Points  PointlList=l  Point2Li8t=2 

Translate  Curve/1  DeltaX=0  DeltaY=0DeltaZ=0, 003526  CID=0  /Points /Create 
Coimect  Curves  CurvelList=l  Curve2List=2  /Points 

Translate  Surface/1  DeltaX=-0 . 01175333  DeltaY=0  DeltaZ=0  CID=0  /Points /Create 
Connect  Surfaces  Surf acelList=2  Surf ace2List— 1  /Points 
MeshParam  Solid/1  Type=Hexa  U=1  V=1  W=3  Pattern-1  PID=1 
Define  Point  X=352.6  Y=-17.63  Z=-17.63  CID=:0 

Translate  Point /9  DelteiX=0  DeltaY=35.26  DeltaZ=0  CID=0  /Create 
Connect  Points  PointlList=9  Point2Li8t=10  . 

Translate  Curve/3  DeltaX=0  DeltaY=0  DeltaZ=35.26  CID=0  /Points/Create 

Connect  Curves  CurvelList=3  Curve2List=4  /Points 

Connect  Surfaces  Surf acelListsl  Surf ace2List=3  /Points 

Me^Param  Solid/2  Type=:Hexa  U=1  V=1  W=996  Pattern:=l  PID=:1  WSpace=500 

Check  Grid/ltlOOOOO  Type=Duplicates  View=l  PartLi0t=O  Tolerance=0 . 000  /  p 

Equivalence  Grids  /Update /NoExtend/NoCollapse  .  .  n  or  3  /Create 

Edit  EOSJWL/1  EOSA=3.712e+ll  EOSB=3 , 231e+09  BOSRl=4.15  EOSR2=0.95  EOSOmega-0 . 3  /Create 

Edit  DetSph/1  MID=1  DetX-0  DetY=0  DetZ=0  DetVel=6390  DetTime=0  i  ■5Q'^7 

Edit  EOSPol/2  EOSAl=2.306e+09  BOSA2=8 . 432e+09  EOSA3=B . 014e+09  EOSB0=0.4934  EOSBl=1.393 

EOSB2=0  EOSB3=Blank  /Create 

Edit  DMat/1  DMatRHO=1630  E0ID=:1  /Create 

Edit  DMat/2  DMatRHO=1025  EOID=2  /Create 

Edit  TICVal/1  Name = "DENSITY"  Value=1630  Pair=l  /Create 

Edit  TICVal/1  Name="SIE"  Value=4290000  Pair=2  /Modify 

Edit  TICVal/2  Name  =  "DENSITY"  Value=:1025  Pair=l  /Create 

Edit  TICVal/2  Name="SIE"  Value=3750.4  Pair=2  /Modify  ^  -1  /Create 

Edit  TICEul/1  Type=Elements  TICEulList=lt3  MID=1  TICValId=l  Level-9  ^ 

Idit  TICEul/l  TWe=Blemeiit8  TICEulLi8t=4t999  MID=2  TICValId=2  Level=l  DMaterial=2  /Modify 
Edit  PEulerl/1  TYPE=MMHydro  TICEulID=l  /Create 

Write  MSCInput  File="bulk.dat"  Format="Bulk"  /NoExec/NoCase/Bulk 


c  VOLUME  TO  RADIUS  CONVERSION  PROGRAM  :  Deep  Spherical  Bubble 

c 

c  program  to  convert  volume  of  a  1-D  wedge  to  spherical  radius 

c 

program  vol2rad 

open(14,file='th_volume_gll.ext ' ,status='old' ) 

open ( 15 ,  f  ile= ' r_vs_t .  ext ' ,  status=: 'new' ) 
do  100  i=l,9 
read(14, *) 

100  continue 

write (*,*)  '#  of  lines  in  volxime  time  history  file?' 
read ( * / * )  n 
n  =  n-9 
do  200  1=1, n 
read (14,*)  t,v 
t  =  t  *  1000 

r  =  ((300.*v)**{l./3.))*100 
write (15,*)  t,r 
200  cont inue 

close (14) 
close (15) 
end 
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c 

c 

c 

c 

c 

c 

c 


200 

100 


400 

300 

999 

998 


INCOMPRESSIBLE /NEGLIGIBLE  INTERNAL  ENERGY  R  VS  T  CALCULATION 


this  program  calculates  r  vs  t  behavior  for  a  spherical  bi^ble 
with  negligible  internal  energy  expanding  and  collapsing  in  an 
incompressible  media 


progreun  void 

implicit  real*16  (a-h,o-z) 

open ( 14, file=' void. ext status unknown  ) 

write {*,*)  "input  r(max):' 
read(*,*)  rm 

write {*,*)  "input  start  time:' 

read(*,*)  tO  o*  \  / 

write(*,*)  'input  #  of  integration  points  (to  max  radius): 

write{*,*)  'input  #  of  integration  steps  per  output  step:' 
read ( * , * )  m 

write{*,*)  'input  water  density:' 
read(*,*)  rho 

write(*,*)  'input  hydrostatic  pressure:' 
read{*,*)  phyd  ^  ^ 

write(*,*)  'input  **tolerance**  factor: 
read{*,*)  tfact 

cinv  =  ((3.*rho)/{2.*phyd))**(0.5) 

k  =  n/m 

n  =  m*k 

rO  =  0.0 

dr  =  (rm-r0)/n 

tol  =  dr/tfact 

write(14,999)  t0»1000, r0*100 
writei*,*) 

period  =  1 . 83*rm* (rho/phyd) ** ( 0 . 5) *1000  ^ 

write (*,998)  'period  from  (8.8)  in  Cole:  ', period 
rm  =  rm+tol 


fO  =  0.0 
do  100  i=l,k 
do  200  j=l,m 

rl  =  r0+( { (i-l) *m)+j)*dr 

fl  =  ((rm/rl)**3.0-1.0)**{-0.5) 

dt  =  ( (fl+f0)/2.0)*dr*cinv 

tl  =  tO+dt 

fO  =  fl 

to  =  tl 

continue 

write(14,999)  tl*1000,rl»100 
cont inue 
do  300  i=l,k 
do  400  j=l,m 

rl  =  rm-(  ( (i-l)*in)+j)  *dr 
£1  =  ( (rm/rl)**3.0-1.0)**(-0.5) 
dt  =  ( (fl+f0)/2.0)*dr*cinv 
tl  =  tO+dt 
fO  =  fl 

to  =  tl 

continue 

write(14,999)  tl*1000,rl*100 
cont inue 
format (2el7 .9) 
format (a, elO . 3 ) 
close ( 14 ) 
end 
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Igure  33.  Bubble  Shape  for  Rigid 


at  h^=4 
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Figure  34.  Bubble  Shape  for  Rigid  Wall  at  h*=4.()15  at  Time  t=0.95T 
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igure  36.  Bubble  Shape  for  Rigid  Wall  at  h*=4.()15  at  Time  t=0.99T 


Figure  52.  Bubble  Shape  for  Rigid  Wall  at  h*-1.000  at  Time  t-0.50T 
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Figure  53.  Bubble  Shape  for  Rigid  Wall  at  h*=  1.000  at  Time  t=0.80T 
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Figure  54.  Bubble  Shape  for  Rigid  Wall  at  h*-1.000  at  Time  t  0.90T 


Figure  55.  Bubble  Shape  for  Rigid  Wall  at  h*=1.000  at  Time  t— 0.95T 
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Figure  56.  Bubble  Shape  for  Rigid  Wall  at  h*=1.000  at  Time  t-0.98T 
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Bubble  Shape  for  Rigid  Wall  at  h*=1.000  at  Time  t=1.00T 


Figure  58. 


APPENDIX  F:  BUBBLES  NEAR  RIGID  BOUNDARIES-INPUT 


$  MSC/DYTRAN  RUN  FILE  :  Bxibble  Near  Rigid  Wall 

$ 

START 
TIME  =999 
CEND 
$ 

TITLE  =  BUBBLE  NEAR  WALL 


$ 

ENDSTEP  =  99999 
ENDTIME  =  17.995E-3 
PARAM,  INISTEP ,  0 . 2E-‘5 
P ARAM , MINSTE P , 0 . 2  E - 1 1 
PARAM, STEPFCT, 19 . 0717446865 
PARAM , DELCLUMP ,0.1 
PARAM , RHOCUT , 1 . E - 4 
PARAM , ROH YDRO , 1 . E - 4 
PARAM ,  ROMULTI ,  1 .  E  --  4 
$ 

TYPE  (MATHIS)  =  TIMEHIS 
MATS  (MATHIS)  =  1 
SET  1=1 

MATOUT  (MATHIS)  =  VOLUME 
STEPS  (MATHIS)  =  0, THRU, END, BY, 1 
SAVE  (MATHIS)  =  99999 
$ 

TYPE  (DENS)  =  ARCHIVE 

ELEMENTS  (DENS)  =  2 

SET  2  =  lt9801 

ELOUT  (DENS)  =  DENSITY 

TIMES  (DENS)  =  0 , THRU, END, BY, 0 . 005 


$ 

BEGIN  BULK 

$ 

TICEUL  1 

+  CYLINDERl 

+  ELEM  2 

CYLINDERl 

+  0.1 

SETl  2  1 

$ 

TICVAL  1 

TICVAL  2 


1 

1 

2 

2 

0.0 

-0.10 

THRU 

23670 

DENSITY 

1630. 

DENSITY 

1025. 

$ 

INCLUDE  bulk.dat 


$ 

ENDDATA 


$ 


9. 

1. 

0.0 


SIE 

SIE 


(Free-Field  Case) 


0.0  0.10  0.0 


4290000. 

20076. 


+ 

+ 


+ 
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MSC/XL  INPUT  FILE  i  Bubble  Near  Rigid  Wall  (Free-Field  Case)  (1/2) 


Edit  ApplicationTable  DuplicateElementTolerance=5e-06 
Edit  ApplicationTable  r>uplicateGridTolerance=5e-06 
Edit  ApplicationTable  DuplicateFEFaceTolerance=5e”06 
Edit  ApplicationTable  DegenerateElementTolerance=5e-06 
Define  Point  X=:0  Y=0  Z-0  CID=0  OutputList=:l 

Translate  Point/1  DeltaX=2.475  DeltaY=0  DeltaZ=0  CID=0  0utputList=2  /Create 
Coimect  Points  PointlList=l  Point2List=2  OutputListsl 

Translate  Curve/1  DeltaX=0  DeltaY=2.475  DeltaZ=0  CID=:0  0utputList=2  /Create/Points 

Connect  Curves  CurvelList=l  Curve2Li8t=2  OutputList=l  /Points 

Define  Point  X=:0.025  Y=0  Z=0  CID=0  0utputList=5 

Define  Point  X=0.05  Y=0  Z=0  CID=0  0utputList=6 

Define  Point  X=:0.075  Y=0  Z=0  CID=0  0utputList=7 

Translate  Point/1,7  DeltaX=0  DeltaY=0.075  DeltaZ=0  CID=0  OutputList=8t9  /Create 

Translate  Point/5t6  DeltaX=0  DeltaY=0.025  DeltaZ=0  CID=0  OutputList=10tll  /Create 

Connect  Points  PointlList=8  Point2List=l  0utputLi8t=3 

Connect  Points  PointlList=10  Point2List=5  0utputList=4 

Connect  Points  PointlList=ll  Point2List=6  0utputList=5 

Connect  Points  PointlList=9  Point2List=7  0utputList=:6 

Connect  Points  PointlList=8  Point2Li8t=9  0utputList=7 

Connect  Points  PointlList=10  Point2List=ll  OutputList=8 

Connect  Curves  CurvelList=3  Curve2List=4  OutputLi8t=2  /Points 

Connect  Curves  CurvelLi8t=4  Curve2List=5  OutputLi8t=3  /Points 

Connect  Curves  CurvelList=8  Curve2Li8t=7  OutputLi8t=4  /Points 

Connect  Curves  CurvelList^B  Curve2Li8t=6  0utputList=5  /Points 

Translate  32  Surface/2t5  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  /Create /NoPoints 
Translate  1  Surf ace/2tl33  DeltaX=0  DeltaY=2.475  DeltaZ=0  CID=0  /Create/NoPoints 
Rotate  Surface/2tl33  FromX=2.475  FromY=0  FromZ=0  ToX=2.475  ToY=0  ToZ=:l  Angler-90 
Of f setAngle=0  CID=0  OutputList=434t441  /Modify/NoPoints  ^ 

Translate  3  Point/4  DeltaX=0.025  DeltaY=0  DeltaZ=0  CID=0  OutputList=12tl4  /Create 

Translate  3  Point /4 , 12tl4  DeltaX=0  DeltaY==0 . 025  DeltaZ:.0  CID=0  OutputLi8t=15t26  /Create 

Connect  Points  PointlList=4  Point2List=23  OutputList=9 

Connect  Points  PointlList=16  Point2List=20  OutputList=10 

Connect  Points  PointlList=17  Point2List=21  OutputList=ll 

Connect  Points  PointlList=14  Point2List=26  OutputList=12 

Connect  Points  PointlList=:23  Point2List=:26  OutputList=13 

Connect  Points  PointlList=20  Point2List=21  OutputList=14 

Connect  Points  PointlList=16  Point2List=17  OutputList=15 

Connect  Points  PointlList=4  Point2List=14  0utputLi8t=16 

Connect  Curves  CurvelList— 9  Curve2List=10  OutputList=266  /Points 

Connect  Curves  CurvelList=13  Curve2List=14  OutputList=267  /Points 

Connect  Curves  CurvelList=14  Curve2List=15  OutputList=268  /Points 

Connect  Curves  CurvelList=15  Curve2List=16  Output List =2 6 9  /Points 

Connect  Curves  CurvelLi8t=ll  Curve2List=12  Output List =2 70  /Points 

Translate  Point /2  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  OutputList=27  /Create 

Translate  Point /3  DeltaX=0  DeltaY=0.075  DeltaZ^O  cp=0  OutputList=28  /Create 

Connect  Points  PointlList=:28  Point2Li8t=26  OutputList=17 

Connect  Points  PointlList=526  Point2List=27  0utputList=18 


Define  Point  X=0  Y=400  Z=0  CID=0  OutputList=29 

Sweep  2  Point /2 9  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=0  ToZ=l  Angle=-45  Of f setAngle-0 


CID=:0 


OutputList=19t20  /Points 

Connect  Curves  CurvelList=17tl8  Curve2List=19t20  OutputList=271t272  /Points 
Translate  Surf ace/lt272  DeltaX=0  DeltaY=0  DeltaZ=-5e-6  CID=0  /Modify/NoPoints 
Translate  Surf ace/lt272  DeltaX==0  DeltaY=0  DeltaZ=10e-6  CID=0  /Create /NoPoints 
Rotate  Surface/lt272  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle-0 


CID=0 


Rotate^Surface/273t544  FromX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=”l  Of f setAngle-0 
CID=0  /Modify/NoPoints 

Connect  Surfaces  Surf acelList=lt272  Surf ace2List=273t544 
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!  MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  (Free -Field  Case)  (2/2) 

I 

MGshParam  Solid/1  Type=Hexa  U=99  V=99  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpacG=l 
ResultingGridListsltlOOOO  OutputLiBt=lt9801  /Extend/Opdate/NoOverwrite/NoMidNode 
Select  Part /I  ^  ^  ^ 

MesbParam  Solid/2t270  Type=Hexa  U=1  V=1  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 
ResultingGridList=20001t22152  OutputList=9802tl0070  /Extend/Update /NoOverwrite/NoMidNode 
MeshParam  Solid/271t272  Type=Hexa  U=34  V=200  W=1  Pattern=l  PID=1  USpacesl  VSpace=120 
WSpace=l  Re8ultingGridList=22153t50292  OutputList=10071t23670 
/ Ext end / Update / NoOve rwr i t e / NoMidNode 

I 

Edit  EOSJWL/1  EOSA=3 .712e+ll  EOSB=3 . 231e+09  EOSRl=4 . 15  EOSR2=0.95  EOSCMega^O . 3  /Create 
Edit  DetSph/1  MID=1  DetX=0  DetY=0  DetZ=0  DetVel=6390  DetTime=0  /Create 

Edit  EOSPol/2  EOSAl=2 .306e+09  EOSA2=8 . 432e+09  EOSA3=8 . 014e+09  EOSB0=0.4934  BOSBl=1.3937 

EOSB2=0  EOSB3=Blank  /Create 

Edit  DMat/1  DMatRHO=1630  EOID=l  /Create 

Edit  DMat/2  DMatRHO=1025  EOID=2  /Create 

Edit  PEulerl/1  TYPE=MMHydro  TICEu 110=1  /Create 

I 

Check  Grid/lt50292  Type=Duplicate8  View=l  PartListsOtl  Tolerance=5e-06  /Update /NoEx tend 
Check  Grid/lt50292  Type=Duplicates  View=l  PartList=Otl  Tolerance=5e-06  /Update /NoEx tend 
Equivalence  Grids  /Update /NoExtend/NoCollapse 

Write  MSCInput  File="bulk . dat "  Formats "Bulk"  /NoExec/NoCase/Bulk 


$ 

$  MSC/DYTRAN  RUN  FILE  ;  Biibble  Near  Rigid  Wall  (h*=4.015) 

$ 

START 
TIME  =  999 
CEND 
$ 

TITLE  =  BUBBLE  NEAR  WALL 


$ 

ENDSTEP  =  9999 
ENDTIME  =  17.995E-3 
PARAM, INISTEP,0.2E-5 
PARAM  r  MINSTEP , 0 . 2E- 11 
PARAM, STEPFCT, 19 . 0717446865 
PARAM, DELCLUMP ,0.1 
PARAM , RHOCUT , 1 . E - 4 
PARAM , ROH YDRO , 1 . E - 4 
P ARAM , ROMULTI , 1 . E - 4 


$ 

TYPE  (MATHIS)  =  TIMEHIS 
MATS  (MATHIS)  =  1 
SET  1=1 

MATOUT  (MATHIS)  =  VOLUME, YMOM 
STEPS  (MATHIS)  =  0, THRU, END, BY, 1 
SAVE  (MATHIS)  =  99999 
$ 

TYPE  (DENS)  =  ARCHIVE 
ELEMENTS  (DENS)  =  2 
SET  2  =  lt21087 
ELOUT  (DENS)  =  DENSITY 
TIMES  (DENS)=  0 007431, . 011890 , 
$ 

BEGIN  BULK 


013376, .014119, .014565, .014713, .014862 


$ 

TICEUL  1 

+  CYLINDERl 

+  ELEM  2 

CYLINDERl 

+  0.1 

SETl  2  1 


119. 

2  2  1. 

0.0  -0.10  0.0  0.0  0.10  0.0 

THRU  42708 


$ 

TICVAL  1 
TICVAL  2 


DENSITY  1630.  SIE  4290000. 

DENSITY  1025.  SIE  20076. 


$ 

INCLUDE  bulk.dat 


$ 

ENDDATA 


$ 
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MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  (h*=4.015)  (1/2) 


Edit  ApplicationTable  DuplicateElemeiitTolerance=5e-06 
Edit  ApplicationTable  DuplicateGridTolerance=5e-06 
Edit  ApplicationTable  I>uplicateFEFaceTolerance=5e-06 
Edit  ApplicationTable  DegenerateElementTolerance^Be-OS 
Define  Point  X=:0  Y=0  Z=0  CID=0  0utputLi8t=l 

Translate  Point/1  DeltaX=2.475  Deltay=0  DeltaZ=0  CID=0  OutputLi8t=:2  /Create 
Connect  Points  PointlLi8t=l  Point2List=2  OutputList=l 

Translate  Curve/1  DeltaX=0  DeltaY=2.475  DeltaZ=0  CID=0  OutputList=2  /Create/Points 

Connect  Curves  CurvelList=l  Curve2List=2  OutputList=:l  /Points 

Define  Point  X=:0.025  Y=0  Z=:0  CID=0  OutputList=5 

Define  Point  X=0.05  Y=0  Z=0  CID=0  OutputList=6 

Define  Point  X=0.075  Y=0  Z=0  CID=0  0utputList=7 

Translate  Point /I, 7  DeltaX=0  DeltaY=0.075  DeltaZ=0  CID=0  OutputList=8t9  /Create 

Translate  Point/5t6  DeltaX=0  DeltaY=0.025  DeltaZ=:0  CID=0  OutputList=10tll  /Create 

Connect  Points  PointlList=:8  Point2List=l  OutputList=3 

Connect  Points  PointlListslO  Point2List=5  OutputLi8t=4 

Connect  Points  PointlList=ll  Point2List=6  0utputList=5 

Connect  Points  PointlList=:9  Point2Li8t=7  0utputList=6 

Connect  Points  PointlList=8  Point2List=9  0utputLi8t=7 

Connect  Points  PointlLi8t=10  Point2List=ll  0utputList=8 

Connect  Curves  CurvelLi8t=3  Curve2Li8t=4  0utputList=2  /Points 

Connect  Curves  CurvelList=:4  Curve2Li8t=5  Output  List =3  /Points 

Connect  Curves  CurvelLi8t=:8  Curve2List=7  0utputList=4  /Points 

Connect  Curves  CurvelList=5  Curve2Li0t=6  0utputList=:5  /Points 

Translate  32  Surface/2t5  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  /Create/NoPoints 
Translate  1  Surf ace/2tl33  DeltaX=0  DeltaY=2.475  DeltaZ=0  CID=0  /Create /NoPoints 
Rotate  Surface/2tl33  FroinX=2,475  FromY=0  FromZ=:0  ToX=2.475  ToY=0  ToZ=l  Angle=-90 
Off8etAngle=0  CID=0  OutputList=434t441  /Modify /NoPoints 

Translate  3  Point/4  DeltaX=0.025  DeltaY=0  DeltaZ=0  CID=:0  OutputList=12tl4  /Create 

Translate  3  Point /4 , 12tl4  DeltaX=0  DeltaY=0.025  DeltaZ=0  CID=0  OutputLi8t=15t26  /Create 

Connect  Points  PointlList=4  Point2List=23  0utputList=9 

Connect  Points  Point lList=:  16  Point2List=20  OutputList=10 

Connect  Points  PointlLi8t=17  Point2List=21  OutputList=ll 

Connect  Points  PointlList=14  Point2List=26  0utputList=12 

Connect  Points  PointlList=23  Point2Li8t=26  0utputList=13 

Connect  Points  PointlLi8t=20  Point2List=21  0utputList=14 

Connect  Points  PointlList=16  Point2List=17  OutputList=15 

Connect  Points  PointlList=4  Point2Li8t=14  0utputList=16 

Connect  Curves  CurvelList-9  Curve2Li8t=10  OutputList=266  /Points 

Connect  Curves  CurvelList=13  Curve2List=14  OutputList=267  /Points 

Connect  Curves  CurvelList=14  Curve2Lists=15  Output  List  =2  6  8  /Points 

Connect  Curves  CurvelList^lB  Curve2Li8t=:16  OutputList=269  /Points 

Connect  Curves  CurvelList=ll  Curve2Li8t=12  Output List =2 70  /Points 

Translate  Point/2  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  OutputList=27  /Create 

Translate  Point /3  DeltaX=0  DeltaY=0.075  DeltaZ=0  CTD=0  OutputList=28  /Create 

Connect  Points  PointlLi8t=28  Point2Li8t=26  OutputList=17 

Connect  Points  PointlLi8t=26  Point2List=27  OutputList=18 


Define  Point  X=0  Y=400  Z=0  CID=0  OutputList=29  i  a 

Sweep  2  Point /29  FromX=0  FroniY=0  FromZ=:0  ToX=0  ToY=0  ToZ=l  Angle=-45  Of f 8etAngle=0 


CIDsO 


OutputList=:19t20  /Points  . 

Connect  Curves  CurvelList=17tl8  Curve2Li8t=19t20  OutputList=271t272  /Points 
Translate  Surf ace/lt272  DeltaX=0  DeltaY=0  DeltaZ=-5e-6  CID=0  /Modify/NoPoints 
Translate  Surf ace/lt272  DeltaX=0  DeltaY=0  DeltaZ=10e-6  CID=0  /Create/NoPoints  ^  ^ 

Rotate  Surface/lt272  FroroX^O  FroinY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle-0  CID-0 


/Modify/NoPoints  ^  ^  ^ 

Rotate  Surface/273t544  FroinX=0  FroinY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=-1  Of f setAngle^O 


CID=0  /Modify/NoPoints 

Connect  Surfaces  Surf acelList=lt272  Surface2List=273t544 
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MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  (h*s:4.015)  (2/2) 

Edit  EOSJWL/1  E0SA=3 . 712e+ll  E0SB=3 . 231e+09  EOSRl=4.15  EOSR2=0.95  EOSOmega=0 . 3  /Create 
Edit  DetSph/1  MID=1  DetX^O  DetY=0  DetZ=0  DetVel=:6390  DetTime=0  /Create 

Edit  EOSPol/2  EOSAl=2.306e+09  EOSA2=8 . 432e+09  EOSA3=8 . 014e+09  EOSB0=0.4934  EOSBl=1.3937 

EOSB2=0  EOSB3=Blank  /Create 

Edit  DMat/1  DMatRHO=1630  EOID=l  /Create 

Edit  DMat/2  DMatRHO=1025  EOID=2  /Create 

Edit  PEulerl/1  TYPE=MMHydro  TICEulID=l  /Create 

t 

Define  Point  X=0  Ys:-2.85  Z=0  CID=0  OutputList=32 

Translate  Point /32  DeltaX=2.475  DeltaY=0  DeltaZ=0  CID=0  OutputList=33  /Create 
Connect  Points  PointlList=32  Point2Li8t=33  OutputList— 21 

Translate  Curve/21  DeltaX=0  DeltaY=2.85  DeltaZ=0  CID=0  OutputLi8t=22  /Points/Create 
Connect  Curves  CurvelList=21  Curve2List=22  OutputList=545  /Points 

Translate  Surface/545  DeltaX=0  DeltaY=0  DeltaZ=-5e-06  CID=0  OutputList =54 5  /Points /Modify 
Translate  Surface/545  DeltaX=0  DeltaY=0  DeltaZ=le-05  CID=0  OutputList =546  /Points/Create 
Rotate  Surface/545  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f 8etAngle=0  CID=0 
OutputList=545  /Points /Modify  *  i  a  n 

Rotate  Surface/546  FromX=0  FroinY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=-1  Of f 8etAngle=0  CID=0 

OutputList =54 6  /Points /Modify  .  .  -  ^  ^  ..t  •  o-m 

Connect  Surfaces  Surf acelList=545  Surf ace2Li8t=546  Weightl=l  Weight2=l  OutputList -2 7 3 


=0  CID=0  OutputList=274t425 


/Points 

Translate  38  Solid/130tl33  DeltaX=0  DeltaY=-0 , 075  DeltaZ= 

/Points/Create 

Define  Point  X=400  Y=0  Z=0  CID=0  OutputList=494 

Translate  Point/494  DeltaX=0  DeltaY=0  DeltaZ=-5e-06  CID=0  OutputList=495  /Create 
Translate  Point/495  DeltaX=0  DeltaY=0  DeltaZ=le-05  CID=0  OutputList=496  /Create 
Rotate  Point/495  FromX=0  FroniY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle=0  CID=0 
Output Li St =4 9 7  /Create  * 

Rotate  Point/496  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=-1  Of f 8etAngle=0  CID=0 

OutputList=498  /Create  . 

Translate  Point/497t498  DeltaX=0  DeltaY=-2.85  DeltaZ=0  CID=0  OutputLi8t=499t500  /Create 
MeshParam  Solid/1  Type=Hexa  0=99  V=99  W=1  Pattern=l  PID=1  USpaee=l  VSpace=l  WSpace=l 
Re8ultingGridList=lt20000  OutputList=lt9801 

MeshParam  Solid/273  Type=Hexa  U=99  V=114  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace-l 
ResultingGridList=20001t43000  OutputList=9802t21087 

MeshParam  Solid/2t270  Type=Hexa  D=1  V=1  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 
Re8ultingGridList=43001t45152  OutputList=21088t21356 

MeshParam  Solid/274t425  Type=Hexa  U=1  V=1  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 
Result ingGridList=45153t463 68  OutputList=21357t21508 

MeshParam  Solid/271t272  Type=Hexa  U=34  V=200  W=1  Pattern=l  PID=1  USpace=l  VSpace=120 
WSpace=l  ResultingGridList=46369t74508  OutputList=21509t35108  ^ 

Connect  Points  PointlLi8t=38  Point2List=41  OutputList=27  /Update /NoOverwrite/Extend 
Connect  Points  PointlList=492  Point2List=493  OutputList =2 8  /Update /NoOverwrite /Extend 
Connect  Points  PointlList=497  Point2Li8t=498  OutputLi8t=29  /Update /NoOverwrite /Extend 
Connect  Points  PointlList=499  Point2List=500  OutputList=30  /Update /NoOverwrite /Extend 
Connect  Curves  CurvelList=27  Curve2List=28  OutputList=549 
/Update /NoOverwrite /Align/Extend/Points 

Connect  Curves  CurvelList=29  Curve2List=30  OutputList=550 

/Update/NoOverwrite/Align/Extend/Points  .  ,  ^  ^  ^  ^ 

Connect  Surfaces  Surf acelList=549  Surf ac62List=550  Weightl=l  Weight2=l  OutputList-426 
/Update/NoOverwrite/Align/Extend/Points 

MeshParam  Solid/426  Type=Hexa  U=1  V=38  W=200  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=120 
Re8ultingGridList=74509t90186  OutputList=35109t42708  /Extend/Update/NoOverwrite/NoMidNode 
Check  Grid/lt90186  Type=Duplicates  View=l  PartList=0tl  Tolerance=5e-06  /Update/NoExtend 
Equivalence  Grids  /Update/NoExtend/NoCollapse  /vt 

Check  Grid/lt90186  Type=Duplicates  View=l  PartList=Otl  Tolerance=5e-06  /Update/NoExtend 
Equivalence  Grids  / Updat e / NoExt end /NoCol lapse 

Write  MSCInput  File='*bulk.dat‘'  Format="Bulk“  /NoExec/NoCase/Bulk 
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$  MSC/DYTRAN  RUN  FILE;  Bubble  Near  Rigid  Wall  (h*=2.008) 

$ 

START 

TIME  =  999 
CEND 
$ 

TITLE  =  BUBBLE  NEAR  WALL 

$ 

ENDSTEP  =  9999 
ENDTIME  17.995E-3 
PARAM, INISTEP, 0 . 2E-5 
PARAM , MINSTEP , 0 . 2E-11 
PARAM, STEPFCT, 19 . 0717446865 
PARAM , DELCLUMP ,0.1 
PARAM , RHOCUT , 1 . E - 4 
PARAM , ROH YDRO , 1 . E - 4 
PARAM , ROMULTI , 1 . E - 4 
$ 

TYPE  (MATHIS)  =  TIMEHIS 
MATS  (MATHIS)  =  1 
SET  1=1 

MATOUT  (MATHIS)  =  VOLUME, YMOM 
STEPS  (MATHIS)  =  0 , THRU, END, BY, 1 
SAVE  (MATHIS)  =  9999 
$ 

TYPE  (DENS)  =  ARCHIVE 
ELEMENTS  (DENS)  =  2 
SET  2  =  ltl5444 
ELOUT  (DENS)  =  DENSITY 

TIMES  (DENS)  =  . 00776 012418 ,. 013971, . 014747 ,. 015057 ,. 015213 , 

$ 

BEGIN  BULK 

$ 

TICEUL  1 


+ 

a. 

CYLINDERl 

1 

1 

9. 

ELEM 

2 

2 

2 

1. 

CYLINDERl 

0,0 

-0.10 

0.0 

0.0 

+ 

0.1 

SETl 

2 

1 

THRU 

33189 

$ 

TICVAL 

1 

DENSITY 

1630. 

SIE 

4290000 

TICVAL 

$ 

2 

DENSITY 

1025. 

SIE 

20076. 

INCLUDE  bulk.dat 
$ 

BNDDATA 

$ 


.015368, .015523 


+ 

+ 


0.0  + 
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I 

!  MSC/XL  INPUT  FILE  ;  Btibble  Near  Rigid  Wall  (h*=2.008)  (1/2) 

I 

Edit  ApplicationTable  DuplicateElementTolerance=5e“06 
Edit  ApplicationTable  DuplicateGridTolerance=Be-06 
Edit  ApplicationTable  DuplicateFEFaceTolerance=5e-06 
Edit  ApplicationTable  DegenerateEleinentTolerance=5e-06 

Define  Point  X=0  y=0  Z=0  CID=0  OutputList=:l  ^  ^ 

Translate  Point/1  DeltaX=2.475  DeltaY=0  DeltaZ=0  CID=0  CmtputList=2  /Create 

Connect  Points  PointlList=l  Point2Li0t=2  OutputLi0t=l  .  o 

Translate  Curve/1  DeltaX=0  DeltaY=2 . 475  DeltaZ=0  CID=0  OutputLr8t=2  /Create/Poxnts 

Connect  Curves  CurvelList=l  Curve2List=2  OutputList=l  /Points 

Define  Point  X=0.025  Y=0  Z=0  CID=0  OutputLi8t=5 

Define  Point  X=0.05  Y=0  Z=0  CID=0  OutputLi8t=6 

Define  Point  X=0.075  Y=0  Z=0  CID=:0  0utputList=7  ^  .  o-q 

Translate  Point/1,7  DeltaX=0  DeltaY=0.075  DeltaZ=0  CID=0  ^ 

Translate  Point/5t6  DeltaX=0  Deltay=0.025  DeltaZ=0  CID=0  OutputList=10tll  /Create 

Connect  Points  PointlList=8  Point2List=l  OutputList— 3 

Connect  Points  PointlList=10  Point2List=5  0utputLi8t=4 

Connect  Points  PointlList=ll  Point2Li8t=6  OutputLi8t=5 

Connect  Points  PointlList=9  Point2List— 7  OutputList— 6 

Connect  Points  PointlList=8  Point2Li0t=9  0utputList=7 

Connect  Points  PointlList=slO  Point2Li8ts:ll  OutputList=8 

Connect  Curves  CurvelList=3  Curve2List=4  0utputList=2  /Points 

Connect  Curves  CurvelList=4  Curve2Li8t=5  0utputList=3  /Points 

Connect  Curves  CurvelList=8  Curve2List=7  OutputList=4  /Points 

Connect  Curves  CurvelList=5  Curve2List=6  0utputList=5  />jr^T>rs^-nt-a 

Translate  32  Surface/2t5  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  /Create/NoPoints 
Translate  1  Surface/2tl33  DeltaX=0  DeltaY=2.475  DeltaZ=0  CID-0 

Rotate  Surface/2tl33  Fron«X=2.475  FromY=0  FromZ=0  ToX=2.475  ToY=0  ToZ=l  Angle— 90 
Off setAngle=0  CID=0  OutputList=434t441  /Modify/NoPoints 

Translate  3  Point /4  DeltaX=0.025  DeltaY=0  DeltaZ=0  CID=0  OutputLi0t=12tl4 

Translate  3  Point/4, 12tl4  DeltaX=:0  DeltaY=0.025  DeltaZ=0  CID=0  OutputList=15t26  /Create 

Connect  Points  PointlList=4  Point2List=23  OutputList-9 

Connect  Points  PointlList=16  Point2Li8t=20  OutputLi8t=10 

Connect  Points  PointlList=17  Point2List=21  Output List =11 

Connect  Points  PointlList=14  Point2List=26  OutputList=12 

Connect  Points  PointlList=23  Point2Li8t=26  0utputList=13 

Connect  Points  PointlList=20  Point2Li8t=21  OutputList=14 

Connect  Points  PointlList=16  Point2List=17  outputList=15 

Connect  Points  PointlList=4  Point2List=14  OutputList-lS 

Connect  Curves  CurvelList=9  Curve2List=10  OutputList=266  /Points 

Connect  Curves  CurvelList=13  Curve2List=14  C)utputList=267  /Points 

Connect  Curves  curvelLi8t=14  Curve2Li8t=15  OutputList=268  /Points 

Connect  Curves  Curvell,i8t=15  Curve2Li8t=16  OutputList=269  /Points 

Connect  Curves  CurvelLi8t=ll  Curve2List=12  OutputLi8t=270  /Points 

Translate  Point/2  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  OutputLi8t=27  /Create 

Translate  Point/3  DeltaX=0  DeltaY=0.075  DeltaZ=0  cro=0  OutputLi8t=28  /Create 

Connect  Points  PointlList=28  Point2List=26  0utputList=17 

Connect  Points  PointlList=26  Point2List=27  0utputList=18 

Define  Point  X=0  Y=:400  Z=0  CID=0  OutputList=29  ,  ,  n  ^  riD-O 

Sweep  2  Point /2 9  FroinX=0  FromY=0  FromZ^O  ToX=0  ToY=0  ToZ=l  Angles -45  OffsetAng 

OutpUtList  =  19t20  /Points  ^  /Points 

Connect  Curves  CurvelList=17tl8  Curve2List=19t20  OutputList=271t272 

Translate  Surf ace/lt272  DeltaX=0  DeltaY=0  DeltaZ=-5e-6  CID=0  /Modify/NoPoints 

Translate  Surfan/It272  DeltaX=0  DeltaY=0  DeltaZ=10e-6  CID=0  /“^“te/NoPoints 

Rotate  Surface/lt272  FromX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle-O  CID-0 

Ro?fte''surface/^3t544  Frona=0  PromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=-1  Of f setAngle=0 


CID=0  /Modify/NoPoints  _ 

Connect  Surfaces  Surf acelLi8t=lt272  Surface2Li8t=273t544 
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!  MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  (h*=2.008)  (2/2) 

I 

Edit  EOSJWL/1  E0SA=3 . 712e+ll  E0SB=3 . 231e+09  E0SR1=4.15  EOSR2=0.95  EOSOmega=0 . 3  /Create 
Edit  DetSph/1  MID=1  DetX=0  DetY=0  DetZ=0  DetVel=6390  DetTime=0  /Create 

Edit  EOSPol/2  E0SA1=2 . 306e+09  EOSA2=8 . 432e+09  EOSA3=8 . 014e+09  EOSB0=0.4934  EOSBl=l,3937 

EOSB2=0  E0SB3=:Blank  /Create 

Edit  DMat/1  DMatRHO=1630  E0ID=1  /Create 

Edit  DMat/2  DMatRHO=1025  E0ID=2  /Create 

Edit  PEulerl/1  TYPE=MMHydro  TICEulID=l  /Create 

I 

Define  Point  X=0  Y=~1.425  Z=0  CID=0  OutputList=32 

Translate  Point/32  DeltaX=2.475  DeltaY=0  DeltaZ=0  CID=0  OutputLi8t=33  /Create 
Connect  Points  PointlList=32  Point2List=33  OutputList=21 

Translate  Curve/21  DeltaX=0  DeltaY=1.425  DeltaZ=0  CID=0  OutputLi8t=22  /Points /Create 
Connect  Curves  CurvelLi8t=21  Curve2List=22  OutputLi8t=545  /Points 

Translate  Surface/545  DeltaX=0  DeltaY^O  DeltaZ=-5e-06  CID=0  OutputList=545  /Points /Modify 
Translate  Surface/545  DeltaX=0  DeltaY=0  DeltaZ=le-05  CID=0  OutputList=546  /Points /Create 
Rotate  Surface/545  PromX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of feet Angle sO  CID=0 
OutputList=545  /Points /Modify 

Rotate  Surface/546  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY^l  ToZ^O  Angle=-1  Of f 8etAngle=0  CID^Q 
OutputList=546  /Points /Modify  . 

Connect  Surfaces  Surf acelLi8t=545  Surf ace2List=546  Weightl=:l  Weight2=l  OutputList=273 
/Points 

Translate  19  Solid/130tl33  DeltaX^O  DeltaY=-0 . 075  DeltaZ-0  CID=0  OutputList=274t349 
/Points/Create 

Define  Point  X=400  Y=0  Z=:0  CID=:0  OutputLi8t=494 

Translate  Point/494  DeltaX=:0  DeltaY=0  DeltaZ=--5e-06  CID=0  Output  List  =4  9  5  /Create 
Translate  Point/495  DeltaX=0  DeltaY^O  DeltaZ=le-05  CID=0  OutputLi8t=496  /Create 
Rotate  Point/495  FromX^O  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f BetAngle=0  CID=0 
OutputList=497  /Create 

Rotate  Point/496  FroinX=0  FromY=0  FromZ:=0  ToX=0  ToY=l  ToZ=0  Angle=-1  Of f 8etAngle=0  CID=0 
Output Li St =498  /Create 

Translate  Point /497t498  DeltaX=0  DeltaY=-l . 425  DeltaZ=0  CID=0  OutputList=499t500  /Create 
MeshParam  Solid/1  Type=Hexa  U=99  V=99  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 
MeshParam  Solid/273  TypG=Hexa  U=99  V=57  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 

Select  Part /I  ^  ^  . 

MeshParam  Solid/2t270  Type=Hexa  U=1  V=1  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 
MeshParam  Solid/274t349  Type=Hexa  U=1  V=1  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 
MeshParam  Solid/271t272  Type=Hexa  U=34  V=200  W=1  Pattern=l  PID=1  USpace=l  VSpace=120 
WSpace=l 

Connect  Points  PointlList=38  Point2List=41  OutputLi8t=27 
Connect  Points  PointlLi8t=264  Point2Li8t=265  OutputList=28 
Connect  Points  PointlList=497  Point2List=498  OutputList=29 
Connect  Points  PointlList=499  Point2Li8t=500  OutputLi8t=30 
Connect  Curves  CurvelLi8t=27  Curve2List=28  OutputList=549  /Points 
Connect  Curves  CurvelLi8t=29  Curve2List=30  Output List =5 50  /Points 

Connect  Surfaces  Surf acelList=549  Surf ace2List=550  Weightl=l  Weight2=l  Output List =42 6 
MeshParam  Solid/426  Type=Hexa  U=1  V=19  W=200  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=pO 
Check  Grid/lt99999  Type=Duplicates  View=l  PartLi8t=0tl  Tolerance=5e-06  / Update /NoExt end 
Equivalence  Grids  / Updat e / NoExt end/ NoCo 11 apse 

Write  MSCInput  File= "bulk.dat"  Format="Bulk"  /NoExec/NoCase/Bulk 
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$  MSC/DYTRAN  RUN  FILE:  Bubble  Near  Rigid  Wall  (h*=1.374) 

$ 

START 

TIME  =  999 
CENT 
$ 

TITLE  :=  BUBBLE  NEAR  WALL 

$ 

ENDSTEP  =  9999 
ENDTIME  =  17.995E-3 
PARAM, INISTEP , 0 . 2E-5 
PARAM, MINSTEP , 0 . 2E-11 
PARAM, STEPFCT, 19 . 0717446865 
PARAM, DELCLUMP ,0.1 
PARAM , RHOCUT , 1 . E - 4 
P ARAM , ROHYDRO , 1 . E - 4 
PARAM , ROMULTI , 1 . E - 4 
$ 

TYPE  (MATHIS)  =  TIMEHIS 
MATS  (MATHIS)  =  1 
SET  1=1 

MATOUT  (MATHIS)  =  VOLUME , YMOM 
STEPS  (MATHIS)  =  0 , THRU, END, BY, 1 
SAVE  (MATHIS)  =  9999 
$ 

TYPE  (DENS)  =  ARCHIVE 
ELEMENTS  (DENS)  =  2 
SET  2  =  ltl3662 


ELOUT 

TIMES 

$ 

BEGIN  : 
$ 

TICEUL 

(DENS) 

(DENS) 

=  DENSITY 
=  .008071, 

.012914, 

.014528, 

.015335, 

.015758, .015819, 

.015981, 

.016142 

BULK 

1 

+ 

+ 

CYLINDERl 

1 

1 

9. 

+ 

4-  ELEM 

CYLINDER! 

2 

2 

0.0 

2 

-0.10 

1. 

0.0 

o 

o 

o 

o 

0.0 

+ 

+ 

0.1 

SETl 

$ 

TICVAL 

2 

1 

THRU 

30183 

1 

DENSITY 

1630. 

SIE 

4290000. 

TICVAL 

2 

DENSITY 

1025. 

SIE 

20076, 

$ 

INCLUDE  bulk.dat 
$ 

ENDDATA 

$ 
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MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  (h*=1.374)  (1/2) 


Edit  ApplicationTable  DuplicateElementTolerance=5e~06 
Edit  ApplicationTable  DuplicateGridTolerance=:5e--06 
Edit  ApplicationTable  DuplicateFEFaceTolerance=5e-06 
Edit  ApplicationTable  DegenerateElementTolerance=5e-06 
Define  Point  X=0  y=0  Z=0  CID=0  OutputLi8t=l 

Translate  Point/1  DeltaX=2.475  DeltaY=0  DeltaZ=0  CID=0  0utputLi8t=2  /Create 
Connect  Points  PointlList=l  Point2Li8t=2  OutputListd 

Translate  Cu3rve/1  Delt6ac=0  DeltaY=2.475  DeltaZ=0  CID=0  OutputLi8t=:2  /Create/Points 

Connect  Curves  CurvelLi8t=l  Curve2List=2  OutputList=l  /Points 

Define  Point  X=0.025  Y=0  Z=0  CID=0  0utputList=5 

Define  Point  X=0.05  Y=0  Z=0  CID=0  0utputList=:6 

Define  Point  X=0.075  Y=0  Z=0  CID=0  0utputList=7 

Translate  Point/1,7  DeltaX=0  DeltaY=0.075  DeltaZ=0  CID=0  OutputLi8t=8t9  /Create 

Translate  Point/5t6  DeltaX=0  DeltaY=0.025  DeltaZ=0  CID=0  Output List =10 til  /Create 

Connect  Points  PointlList=8  Point2Li8t=l  OutputList=3 

Connect  Points  PointlList=10  Point2List=5  OutputList=4 

Connect  Points  PointlList=ll  Point2List=6  OutputList=5 

Connect  Points  PointlLi8t=9  Point2List=7  OutputLi8t=6 

Connect  Points  PointlLi8t=e  Point2Li8t=9  OutputLi8t=7 

Connect  Points  PointlList=10  Foint2List=ll  0utputList=8 

Connect  Curves  CurvelList=3  Curve2Li8t=4  Output Lis t=2  /Points 

Connect  Curves  CurvelLi8t=4  Curve2List=5  Output List =3  /Points 

Connect  Curves  CurvelLi8t=8  Curve2List=7  OutputList=4  /Points 

Connect  Curves  CurvelList=5  Curve2List=6  OutputList=5  /Points 

Translate  32  Surface/2t5  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  /Create /NoPoints 
Translate  1  Surf ace/2tl33  DeltaX=0  DeltaY=2.475  DeltaZ=0  CID=0  /Create /NoPoints 
Rotate  Surface/2tl33  FromX=2.475  Fromy=0  FromZ=0  ToX=2.475  ToY=0  ToZ=l  Angle=-90 
OffsetAngle=0  CID=0  OutputList=434t441  /Modify/NoPoints 

Translate  3  Point /4  DeltaX=0.025  DeltaY=0  DeltaZ=0  CID=0  OutputLi8t=12tl4  /Create 

Translate  3  Point/4 , 12tl4  DeltaX=0  DeltaY=0.025  DeltaZ=0  CID=0  OutputList=15t26  /Create 

Connect  Points  PointlList=4  Point2Li8t=23  OutputLi8t=9 

Connect  Points  PointlList=16  Point2List=20  OutputListslO 

Connect  Points  PointlList=17  Point2Li8t=21  OutputList=ll 

Connect  Points  PointlList=14  Point2List=26  0utputList=12 

Connect  Points  PointlList=23  Point2Li8t=26  OutputList=13 

Connect  Points  PointlLi8t=20  Point2Li8t=21  OutputList=14 

Connect  Points  PointlList=16  Point2List=17  0utputList=15 

Connect  Points  PointlList=4  Point2List=14  0utputLi8t=16 

Connect  Curves  CurvelList=9  Curve2List=10  OutputLi8t=266  /Points 

Connect  Curves  CurvelLi8t=13  Curve2List=14  OutputList=267  /Points 

Connect  Curves  CurvelLi8t=14  Curve2List=15  OutputList=268  /Points 

Connect  Curves  CurvelLi8t=15  Curve2Li8t=16  OutputLi8t=269  /Points 

Connect  Curves  CurvelList=ll  Curve2Li8t=12  OutputList=270  /Points 

Trzinslate  Point /2  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  OutputLi8t=27  /Create 

Translate  Point/3  DeltaX=0  DeltaY=0.075  DeltaZ=0  CID=0  OutputList=28  /Create 

Connect  Points  PointlLi8t=28  Point2List=26  OutputLi8t=17 

Connect  Points  PointlList=26  Point2Li8t=27  OutputList=18 

Define  Point  X=0  Y=400  Z=0  CID=0  OutputList=29 

Sweep  2  Point/29  FromX=0  FromY=0  FromZ=0  ToX=0  ToY=0  ToZ=l  Angle=-45  Of f 8etAngle=0  CID=0 
OutputList=19t20  /Points 

Connect  Curves  CurvelList=17tl8  Curve2List=19t20  OutputList=271t272  /Points 

Tremslate  Surf ace/lt272  DeltaX=0  DeltaY=0  DeltaZ=-5e-6  CID=0  /Modify/NoPoints 

Translate  Surf ace/lt272  DeltaX=0  DeltaY=0  DeltaZ=10e-6  CID=0  /Create /NoPoints 

Rotate  Surface/lt272  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle=0  CID=0 

/Modify/NoPoints 

Rotate  Surface/273t544  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=’l  Of f setAngle=0 
CID=0  /Modify/NoPoints 

Connect  Surfaces  Surf acelList=lt272  Surf ace2List=273t544 
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I 

!  MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  (h*=1.374)  (2/2) 

Edit  EOSJWL/1  EOSA=3.712e+ll  E0SB=3 . 231e+09  EOSRl=4.15  EOSR2=0.95  EOSOmega=0 . 3  /Create 
Edit  DetSph/1  MID=1  DetX=0  DetY=0  DetZ=:0  DetVel=6390  DetTimG  =  0  /Create 

Edit  EOSPol/2  E0SA1=2 .306e+09  EOSA2=8 . 432e+09  EOSA3=8 . 014e+09  EOSB0=0.4934  EOSBl=1.3937 

EOSB2=0  E0SB3=BlaBk  /Create 

Edit  DMat/1  DMatRHO=1630  E0ID=1  /Create 

Edit  DMat/2  DMatRHO=:1025  E0ID=2  /Create 

Edit  PEulerl/l  TYPE=MMHydro  TICEulID=l  /Create 

I 

Define  Point  X=0  Y=-0.975  Z=0  CID=0  OutputList=32 

Translate  Point/32  DeltaX=2.475  DeltaY=0  DeltaZ=0  CID=0  OutputList=33  /Create 
Connect  Points  PointlList=32  Point2List=33  OutputList=21  ,  .  ^ 

Translate  Curve/21  DeltaX=0  DeltaY=0.975  DeltaZ=0  CID=0  OutputList=22  /Points/Create 
Connect  Curves  CurvelList=21  Curve2List=22  OutputList=545  /Points  fv 

Translate  Surface/545  DeltaX=0  DeltaY=0  DeltaZ=-5e-06  CID=0  0"tputLi8t=545  /Points/M^lfy 
Translate  Surface/545  DeltaX=0  DeltaY=0  DeltaZ=le-05  CID=0  ^tputLlst=546  /Points/Create 
Rotate  Surface/545  ProinX=0  FromY=0  PromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle-0  CID- 

OutputList=545  /Points /Modify  *  r'rn-n 

Rotate  Surface/546  FromX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angler-1  Of f setAngle-O  CID-0 

OutputList=546  /Points /Modify  ,, 

Connect  Surfaces  Surf acelList=545  Surface2List=546  Weightl=l  Weight2=l  OutputLlst-273 

Translate  13  Solia/130tl33  DeltaX=0  DeltaY=-0 . 075  DeltaZ=0  CID=0  0utputList=274t349 
/Points /Create 

Define  Point  X=400  Y=0  Z=0  CID=0  OutputList=494 

Translate  Point/494  DeltaX=0  DeltaY=0  DeltaZ=-5e-06  CID=0  OutputList=495  /Create 
Translate  Point/495  DeltaX=0  DeltaYrO  DeltaZ=le-05  CID=0  ^tputLi8t=496  /Create 
Rotate  Point /495  FromX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle-0  CID-0 

Rot?te^Point%96*^Pro^iX=0  FroinY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=-1  Of f 8etAngle=0  CID=0 

Translate*^PoiL/497t498  DeltaX=0  DeltaY=-l  .425  DeltaZ=0  CID=0  OutputLi8t=499t500  /Create 

MeshParam  Solid/1  Type=Hexa  0=99  V=99  W=1  Pattern=l  PID=1  ySpace=l  ySpace=l 

MeshParam  Solid/273  Type=Hexa  0=99  V=39  W=1  Pattern=l  PID=1  0Space=l  VSpace-1  WSpace-1 

Me8ML^''soiid/2t270  Type=Hexa  0=1  V=1  W=1  Pattern=l  PID=1  0Space=l  VSpace=l 
MeshParam  Solid/274t349  Type=Hexa  0=1  V=1  W=1  Pattern=l 

MeshParam  Solid/271t272  Type=Hexa  0=34  V=200  W=1  Pattern=l  PID=1  0Space=l  VSpace-120 


WSpace=l 

Connect  Points  PointlList=38  Point2List=41  OutputList=27 
Connect  Points  PointlList=192  Point2List=193  OutputList=28 
Connect  Points  PointlLi8t=497  Point2List=498  OutputList=29 
Connect  Points  PointlList=499  Point2Li8t=500  OutputLi8t=30 
Connect  Curves  CuirvelList=27  Curve2List=28  OutputLi8t=549  /Points 
Connect  Curves  CurvelList=29  Curve2LiBt=30  OutputLi8t:=550  /Points 
Connect  Surfaces  Surf acelLi8t=549  Surf ace2List=550  Weightl=l  Weight2-1 
MeshParam  Solid/426  Type=Hexa  U=1  V=13  W=200  Pattern=l  PID=1  USpace=l 
Check  Grid/lt99999  Type^Duplicates  View=l  PartList=Otl  Tolerance=5e-06 
Equivalence  Grids  /Update/NoExtend/NoCollapse 

Write  MSCInput  Filess"bulk.dat"  Formats: "Bulk"  /NoExec/NoCase/Bulk 


Output List s42 6 
VSpace=l  WSpace=120 
/Update /NoExt end 
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$  MSC/DYTRAN  RUN  FILE:  Bubble  Near  Rigid  Wall  (h*=1.000> 

$ 

START 
TIME  =  999 
CEND 
$ 

TITLE  =  BUBBLE  NEAR  WALL 


$ 

ENDSTEP  =  9999 
ENDTIME  =  17.995E-’3 
PARAM, INISTEP, 0 . 2E-5 
PARAM , MINSTEP , 0 . 2 E - 1 1 
PARAM, STEPFCT, 19 . 0717446865 
PARAM , DELCLUMP ,0.1 
PARAM ,  RHOCUT ,  1 .  E  -  4 
PARAM, ROHYDRO, 1 . E~4 
PARAM, ROMULTI , 1 . E-4 
$ 

TYPE  (MATHIS)  =  TIMEHIS 
MATS  (MATHIS)  =  1 
SET  1=1 

MATOUT  (MATHIS)  =  VOLUME, YMOM 
STEPS  (MATHIS)  =  0 , THRU, END, BY, 1 
SAVE  (MATHIS)  =  9999 
$ 

TYPE  (DENS)  =  ARCHIVE 
ELEMENTS  (DENS)  =  2 
SET  2  =  ltl9602 
ELOUT  (DENS)  =  DENSITY 

TIMES  (DENS)  =. 008522 ,. 013636 ,. 015340 , 

$ 

SPC=1 


$ 

BEGIN  BULK 

$ 

TICEUL  1 

+  CYLINDERl 

-I-  ELEM  2 

CYLINDERl 

+  0.1 

SETl  2  1 

$ 

TICVAL  1 

TICVAL  2 

$ 

INCLUDE  bulk.dat 
INCLUDE  wall.dat 
$ 

ENDDATA 

$ 


1 

1 

2 

2 

0.0 

o 

O 

THRU 

47340 

DENSITY 

1630. 

DENSITY 

1025. 

016193, ,016704, .016875, 


9. 

1. 

0.0  0.0  0.10 


SIE  4290000. 

SIE  20076. 


.017045 


+ 

+ 

0.0  + 
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$  MSC/DYTRAN  DATA  FILE  "wall.dat":  Bxibble  Near  Rigid  Wall  (h*  =  1.000) 
$ 

BEGIN  BULK 
$ 


$  THIS 
$ 

SECTION 

CONTAINS 

BULK  DATA 

$ 

GRID 

100585 

-.025 

-.7098 

-200. 

GRID 

100586 

400.025 

-.7098 

-200. 

GRID 

100587 

-.025 

-.7098 

200. 

GRID 

100588 

400.025 

-.7098 

200. 

GRID 

100589 

-.025 

-400.71 

-200. 

GRID 

100590 

400.025 

-400.71 

-200. 

GRID 

100591 

-.025 

-400.71 

200. 

GRID 

100592 

400.025 

-400.71 

200. 

$ 

CHEXA 

47341 

3 

100585 

100589 

100590 

+ 

$ 

$  THIS 
$ 

100592 

100588 

SECTION 

CONTAINS 

THE  LOADS ,  CONSTRAINTS 

$ 

SPC 

1 

100585 

123456 

SPC 

1 

100586 

123456 

SPC 

1 

100587 

123456 

SPC 

1 

100588 

123456 

SPC 

1 

100589 

123456 

SPC 

1 

100590 

123456 

SPC 

1 

100591 

123456 

SPC 

1 

100592 

123456 

100586  100587  100591  + 


THIS  SECTION  CONTAINS  THE  DEFINED  FEFACES  OF  ELEMENTS 


$ 

$ 

$ 

$ 

CFACE 

CFACE 

CFACE 

CFACE 

CFACE 

CFACE 

$ 

PSOLID 

$ 

DMATEL 

$ 

SURFACE 

$ 

COUPLE 

ENDDATA 


7890. 


47341 

47341 

47341 

47341 

47341 

47341 


1.95+11  .28 
SEG  2 

Inside  On 


On 


no 


!  MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  {h*=1.000)  (1/2) 

I 

Edit  ApplicationTable  DuplicateBleinentToleraiice=5e-06 
Edit  ApplicationTable  DuplicateGridTolerance=5e-06 
Edit  ApplicationTable  DuplicateFEPaceTolerance=5e-06 
Edit  ApplicationTable  DegenerateElementTolerance=5e-06 
Define  Point  X=0  Y=0  Z=0  CID=0  0utputLi8t=l 

Translate  Point /I  DeltaX=2.475  Deltay=0  DeltaZ=0  CID=0  0utputList=2  /Create 
Connect  Points  PointlList=l  Point2List=2  0utputLi8t=l 

Translate  Curve/1  DeltaX=0  DeltaY=:2 . 475  DeltaZ^O  CID=0  OutputLi8t=2  /Create/Points 

Connect  Curves  CurvelList=:l  Curve2Li8t=2  OutputList^l  /Points 

Define  Point  X=0.025  Y=0  Z=0  CID=0  0utputList=5 

Define  Point  X=0.05  Y=0  Z=:0  CID=0  OutputLi8t=6 

Define  Point  X=:0.075  Y=0  Z=0  CID=0  0utputList=7 

Translate  Point /I, 7  BeltaX=0  DeltaYsO . 075  DeltaZ=0  CID=0  OutputLi8t=8t9  /Create 

Translate  Point /5t 6  DeltaX=0  DeltaY=0.025  DeltaZ=0  CID=0  OutputLi8t=10tll  /Create 

Connect  Points  PointlList=8  Point2Li8t=l  Output List =3 

Connect  Points  PointlList=10  Point2List=5  OutputList=4 

Connect  Points  PointlLi8t=ll  Point2Li8t=6  OutputList=5 

Connect  Points  PointlLi8t=9  Point2List=7  OutputLi8t=6 

Connect  Points  PointlList=8  Point2List=9  OutputList=7 

Connect  Points  PointlList=10  Point2List=ll  OutputList=8 

Connect  Curves  CurvelList=3  Curve2List=4  OutputList=:2  /Points 

Connect  Curves  CurvelList=4  Curve2Li8t=5  OutputList=3  /Points 

Connect  Curves  CurvelLi8t=:8  Curve2List=:7  OutputLi8t=4  /Points 

Connect  Curves  CurvelList=5  Curve2List=6  OutputList=5  /Points 

Translate  32  Surface/2t5  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  /Create/NoPoints 
Translate  1  Surf ace/2tl33  DeltaX=0  DeltaY=2.475  DeltaZ=0  CID-0  /Create/NoPoints 
Rotate  Surface/2tl33  FromX=2.475  FroinY=0  FromZ=:0  ToX=2.475  Toy=0  ToZ=l  Angle=-90 
Of fsetAngle^O  CID=0  OutputList=434t441  /Modify/NoPoints 

Translate  3  Point /4  DeltaX=0.025  DeltaY=0  DeltaZ=0  CID=0  Output List =12 t 14  /Create 

Translate  3  Point/4, 12tl4  DeltaX=0  DeltaY=0.025  DeltaZ=0  CID=0  OutputLi8t=15t26  /Create 

Connect  Points  PointlList=4  Point2List=23  OutputList=9 

Connect  Points  PointlList=16  Point2List=20  OutputList=10 

Connect  Points  PointlLi8t=17  Point2List=21  0utputLi8t=ll 

Connect  Points  PointlList=14  Point2List=26  0utputList=12 

Connect  Points  PointlLi8t=23  Point2Li8t=26  OutputList=13 

Connect  Points  PointlLi8t=20  Point2Li8t=21  OutputList=14 

Connect  Points  PointlList=16  Point2Li8t=17  OutputList=15 

Connect  Points  PointlLi8t=4  Point2LiBt=14  OutputList=16 

Connect  Curves  CurvelLi8t=9  Curve2LiBt=10  OutputLi8t=266  /Points 

Connect  Curves  CurvelList=13  Curve2List=14  Output List =2 6 7  /Points 

Connect  Curves  CurvelLi8t=14  Curve2List=15  OutputList=268  /Points 

Connect  Curves  CurvelList=15  Curve2List=16  OutputLi8t=269  /Points 

Connect  Curves  CurvelList=ll  Curve2Li8t=12  OutputList=270  /Points 

Translate  Point/2  DeltaX=0.075  DeltaY=0  DeltaZ=0  CID=0  OutputList=27  /Create 

Translate  Point/3  DeltaX=0  DeltaY=0.075  DeltaZ=0  CID=0  OutputList=28  /Create 

Connect  Points  PointlLi8t=28  Point2Li8t=26  0utputList=17 

Connect  Points  PointlList=26  Point2Li8t=27  OutputList=18 

Define  Point  X=0  Y=400  Z=0  CID=0  OutputList=29 

Sweep  2  Point/29  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=0  ToZ=l  Angle=-45  Of f setAngle=0  CID=0 
OutputLi8t=19t20  /Points 

Connect  Curves  CurvelList=17tl8  Curve2Li8t=19t20  0utputLi8t=271t272  /Points 

Translate  Surf ace/lt272  DeltaX=0  DeltaY=0  DeltaZ=-5e-6  CID=0  /Modify/NoPoints 

Translate  Surface/lt272  DeltaX=0  DeltaY=0  DeltaZ=10e-6  CID=0  /Create/NoPoints 

Rotate  Surf ace/lt272  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=l  Of f setAngle=0  CID=0 

/Modify/NoPoints 

Rotate  Surface/273t544  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l  ToZ=0  Angle=-1  Of f setAngle=0 
CID=0  /Modify/NoPoints 

Connect  Surfaces  Surf acelLi8t=lt272  Surf ace2List=273t544 
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!  MSC/XL  INPUT  FILE  :  Bubble  Near  Rigid  Wall  (h*=1.000)  (2/2) 

I 

Ref lect/NoOverwrite/Extend/Update/Create  Solid/lt272  FroinX=0  FromY=0  FromZ=0  ToX=0  ToY=l 
ToZ=0  CID=0  OutputList=273t544  /NoPoints/NoFlip 

MeshParam  Solid/1,273  Type=Hexa  U=99  V=99  W=1  Pattern=l  PID=1  USpace=l  VSpace=l  WSpace=l 
Result iiigGridList=lt40 000  Output Lis t=lt 19 602  /Extend/Update /NoOverwrite/NoMidNode 
Select  Part /I 

MeshParam  Solid/2 t27 0, 274t 542  Type=Hexa  U=1  V=:l  W=1  Pattern=l  PID=1  USpace=l  VSpace=l 
WSpace=l  Result ingOridList=:40001t443 04  OutputList=19603t20140 
/Extend/Update /NoOverwr it e /NoMidNode 

MeshParam  Solid/271,272,543,544  Type=Hexa  U=34  V=200  W=1  Patterii=l  PID=1  USpace=l 
VSpace=120  WSpace=l  ResultingGridList=44305tl00584  OutputList=20141t47340 
/Extend /Update/ NoOve  rwr i t  e / NoMidNode 

Check  Grid/ltl00584  Type =I>upli cates  View=l  PartList^Otl  Tolerance=5e-06  /Update /NoExt end 
Equivalence  Grids  /Update /NoExt end/NoColl apse 

Check  Grid/ It 100583  Type=Duplicates  View=l  PartList=0tl  Tolerance=5e-06  /Update /NoExt end 
Equivalence  Grids  /Update /NoExt end/NoCol lapse 

I 

Edit  EOSJWL/1  EOSA=3 .712e+ll  EOSB=:3 . 231e+09  EOSRl=:4.15  EOSR2=0.95  EOSOmega=0.3  /Create 
Edit  DetSph/1  MID=1  DetX=:0  DetY=0  DetZ=0  DetVel=6390  DetTime=0  /Create 

Edit  EOSPol/2  EOSAl=2 . 306e+09  EOSA2=8 . 432e+09  EOSA3=8 ♦ 014e+09  EOSB0=0.4934  EOSBl=1.3937 

EOSB2=0  EOSB3=Blank  /Create 

Edit  DMat/1  DMatRHO=1630  EOID-1  /Create 

Edit  DMat/2  DMatRHO=1025  EOID=:2  /Create 

Edit  PEulerl/1  TYPE=MMHydro  TICEulID=l  /Create 

Write  MSCInput  File="bulk.dat "  Format="Bulk"  /NoExec/NoCase/Bulk 


112 


c  VOLUME  TO  RADIUS  CONVERSION  PROGRAM  :  Bubble  Near  Rigid  Wall 

c  { ax i symmetric  1/2  model) 

c 

program  vol2rad 

open  ( 14 ,  files '  th_volume__gll .  ext '  ,  statuss 'old' ) 
open(15,file='r_vs_t.ext' ,status='new' ) 

do  100  i=l,9 
read(14,*) 

100  continue 

write(*,*)  '#  of  lines  in  volume  time  history  file?' 
read {* ,*)  n 
n  =  n-'9 
do  200  i=l,n 
read (14,*)  t,v 
t  =  t  *  1000 

r  =  ((0.5*85.94367*v)**(l./3.))*100 
write (15,*)  t,r 
200  continue 

close (14) 
close (15) 
end 
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c 

c 

c 


C.M.  DISPLACEMENT  CALCULATION  PROGRAM:  Bubble  Near  Rigid  Wall 


progreim  ydisp 

implicit  double  precision (a-h, o-z) 
open {14, file='th_ymom_gll.ext " , status^ "old' ) 
open (15, f ile= "ydisp_vs_t .ext ' , status='new" ) 
do  100  i=l,9 
read{14,*) 

100  continue 

write (*,*)  "#  of  lines  in  momentum  time  history  file?" 

read (* ,*)  n 
n  =  n-9 

w  =  10.2415920507/180.0 
sum  =  0.0 

read (14,*)  tl,ymoml 
write (15,*)  sum, sum 
do  200  i=l,n-l 
read (14,*)  t2,ymom2 

sum  =  sum+(t2-tl)  *0.5*  {ymom2+ymoml) 
t  =  tl+0.5*(t2-tl) 
y  =  sum/w 

write(15,*)  t*1000,y*100 
tl  =  t2 
ymoml  =  ymom2 
200  continue 
close (14) 
close (15) 
end 
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